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ABSTRACT 

We develop a general framework for interpreting and analyzing high-precision 
lightcurves from eccentric stellar binaries. Although our methods are general, we 
focus on the recently discovered Kepler system KOI-54, a face-on binary of two A 
stars with e = 0.83 and an orbital period of 42 days. KOI-54 exhibits strong el- 
lipsoidal variability during its periastron passage; its lightcurve also contains ~ 20 
pulsations at perfect harmonics of the orbital frequency, and another ~ 10 nonhar- 
monic pulsations. Analysis of such data is a new form of asteroseismology in which 
oscillation amplitudes and phases rather than frequencies contain information that 
can be mined to constrain stellar properties. We qualitatively explain the physics of 
mode excitation and the range of harmonics expected to be observed. To quantita- 
tively model observed pulsation spectra, we develop and apply a linear, tidally forced, 
nonadiabatic stellar oscillation formalism including the Coriolis force. We produce 
temporal power spectra for KOI-54 that are semi-quantitatively consistent with the 
observations. Both stars in the KOI-54 system are expected to be rotating pseudosyn- 
chronously, with resonant nonaxisymmetric modes providing a key contribution to 
the total torque; such resonances present a possible explanation for the two largest- 
amplitude harmonic pulsations observed in KOI-54, although we find problems with 
this interpretation. We show in detail that the nonharmonic pulsations observed in 
KOI-54 can be explained by nonlinear thrcc-mode coupling. The methods developed 
in this paper can be generalized in the future to determine the best-fit stellar param- 
eters given pulsation data. We also derive an analytic model of KOI-54's ellipsoidal 
variability, including both tidal distortion and stellar irradiation, which can be used 
to model other similar systems. 

Key words: binaries: close - asteroseismology - stars: oscillations - hydrodynamics 
- waves 



1 INTRODUCTION 



The recently discovered Kepler system KOI-54 ( |Welsh et al. 



2011 henceforth Wll) is a highly eccentric stellar binary 
with a striking lightcurve: a 20-hour 0.6% brightening oc- 
curs with a periodicity of 41.8 days, with lower-amplitude 
perfectly sinusoidal oscillations occurring in between. Such 
observations were only possible due to the unprecedented 
photometric precision afforded by Kepler. Wll arrived at 
the following interpretation of these phenomena: during the 
periastron passage of the binary, each of its two similar A 
stars is maximally subjected to both its companion's tidal 
force and radiation field. The tidal force causes a prolate 
ellipsoidal distortion of each star known as the equilibrium 
tide, so that the resulting perturbations to both the stellar 
cross section and the emitted stellar flux produce a change 



in the observed flux. Along with the effects of irradiation, 
this then creates the large brightening at periastron. Sec- 
ondly, the strong tidal force also resonantly excites stellar 
eigenmodes during periastron, which continue to oscillate 
throughout the binary's orbit due to their long damping 
times; this resonant response is known as the dynamical 
tide. 

Wll successfully exploited KOI-54's periastron flux 
variations, known traditionally as ellipsoidal variability, by 
optimizing a detailed model against this component of 
KOI-54's lightcurve ( |Orosz fc Hauschildt|2000| . In this way, 
Wll were able to produce much tighter constraints on stellar 
and orbital parameters than could be inferred through tradi- 
tional spectroscopic methods alone. Wll also provided data 
on the dynamical tide oscillations. Thirty such pulsations 
were reported, of which roughly two-thirds have frequencies 
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at exact harmonics of the orbital frequency. It is the analy- 
sis of these and similar future data that forms the basis of 
our work. 

In close binary systems, tides provide a key mechanism 
to circularize orbits and synchronize stellar rotation with 
orbital motion. An extensive literature exists on the the- 
ory of stellar tides (e.g., Zahn 1975 Goodman & Dickson 
|1998| |Witte fc Savonije|[l999j . We have synthesized this 
theoretical formalism, together with other aspects of stellar 
oscillation theory, in order to model the dynamical tide of 
KOI-54 as well as to provide a framework for interpreting 
other similar systems. 

The methods we have begun to develop are a new form 
of asteroseismology, a long-standing subject with broad util- 
ity. In traditional asteroseismology, we observe stars in 
which internal stellar processes (e.g., turbulent convection 
or the kappa mechanism) drive stellar eigenmodes, allowing 
them to achieve large amplitudes ( [Christensen-Dalsgaard| 
|2003[ ). In this scenario, modes ring at their natural frequen- 
cies irrespective of the excitation mechanism. The observed 
frequencies (and linewidths) thus constitute the key infor- 
mation in traditional asteroseismology, and an extensive set 
of theoretical techniques exist to invert such data in order to 
infer stellar parameters and probe different aspects of stellar 
structure ( Unno et al.|1989| . 

In tidal asteroseismology of systems like KOI-54, how- 
ever, we observe modes excited by a periodic tidal poten- 
tial from an eccentric orbit; tidal excitation occurs predomi- 
nantly at I = 2 (§ 3.2 1. Since orbital periods are well below a 
star's dynamical timescale, it is g-modes (buoyancy waves) 
rather than higher-frequency p-modes (sound waves) that 
primarily concern us. Furthermore, since modes in our case 
are forced oscillators, they do not ring at their natural eigen- 
frequencies, but instead at pure harmonics of the orbital fre- 
quency. (We discuss nonharmonic pulsations in § 6.5 ) It is 



thus pulsation amplitudes and phases that provide the key 
data in tidal asteroseismology. 

This set of harmonic amplitudes and phases in principle 
contains a large amount of information. One of the goals for 
future study is to determine exactly how the amplitudes can 
be optimally used to constrain stellar properties, e.g., the 
radial profile of the Brunt- Vaisala frequency. In this work, 
however, we focus on the more modest tasks of delineating 
the physical mechanisms at work in eccentric binaries and 
constructing a coherent theoretical model and corresponding 
numerical method capable of quantitatively modeling their 
dynamical tidal pulsations. 

This paper is organized as follows. In § [2] we give es- 
sential background on KOI-54. In §[3] we give various the- 
oretical results that we rely on in later sections, including 
background on tidal excitation of stellar eigenmodes (§ |3.2[ ), 
techniques for computing disk-averaged observed flux per- 



turbations (§ 3.3 1, and background on including the Coriolis 
force using the traditional approximation (§ 3.4|. In §[4] we 
use these results to qualitatively explain the pulsation spec- 
tra of eccentric stellar binaries, particularly what governs 
the range of harmonics excited. In § [5] we confront the ro- 
tational evolution of KOI-54's stars, showing that they are 
expected to have achieved a state of stochastic pseudosyn- 
chronization. 

In § [6] we present the results of our more detailed mod- 
eling. This includes an analytic model of ellipsoidal variabil- 



Table 1. List of KOI-54 system parameters as determined by 
Wll. Selected from Table 2 of Wll. The top rows contain stan- 
dard observables from stellar spectroscopy, whereas the bottom 
rows result from Wll's modeling of photometric and radial ve- 
locity data. Symbols either have their conventional definitions, 
or are defined in § |3.1| Note that Wll's convention is to use the 
less-massive star as the primary. 
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ity (§ 6.1 1, the effects of nonadiabaticity (§ 6.2 1, the effects 



of fast rotation (§ 
01 



6.3l 



and a preliminary optimization of 
our nonadiabatic method against KOI-54's pulsation data 
We show in 



6.4 1 



6.5 



that the observed nonharmonic 
pulsations in KOI-54 are well explained by nonlinear three- 
mode coupling, and perform estimates of instability thresh- 
olds, which may limit the amplitudes modes can attain. We 
also address whether the highest-amplitude observed har- 
monics in KOI-54 are signatures of resonant synchronization 
locks in § |6.6| We present our conclusions and prospects for 
future work in § [7] 

A few weeks prior to the completion of this manuscript, 
we became aware of a complementary study of KOI-54's 



pulsations (Fuller & Lai||2011 1 



2 BACKGROUND ON KOI-54 

Table [1] gives various parameters for KOI-54 resulting from 
Wll's observations and modeling efforts. Table[2]gives a list 
of the pulsations Wll reported, including both frequencies 
and amplitudes. 



2.1 Initial rotation 

KOI-54's two components are inferred to be A stars. Iso- 
lated A stars are observed to rotate much more rapidly than 
e.g. the Sun, with typical surface velocities of ~ 100 km/s 
and rotation periods of ~ 1 day (Adclman 2004). This re- 
sults from their lack of a significant convective envelope, 
which means they experience less-significant magnetic brak- 
ing, allowing them to retain more of their initial angular 
momentum as they evolve onto the main sequence. We thus 
operate under the assumption that both component stars of 
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Table 2. Thirty largest KOI-54 pulsations. Originally Table 3 
from Wll. Asterisks (*) denote pulsations which are not obvious 
harmonics of the orbital frequency. 
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KOI-54 were born with rotation periods of roughly 
-Pbirth ~ 1 day. 

2.2 Rotational inclination 

Wll constrained both stars' rotation, via line broadening, 
to be f r ot sin i* = 7.5 ± 4.5 km/s (the same for both stars). 
Using the mean values of Ri and R2 obtained from Wll's 
modeling, we can translate this into the following constraints 
on rotation periods (in days): 



9.2 < Pi/sinii < 37 



and 



9.8 < P 2 /sini 2 < 39, 



where (11,12) and (Pi,p2) are the rotational angular mo- 
mentum inclinations with respect to the observer and the 
stellar rotation periods, respectively. If we assume that tidal 
interactions cause both stellar rotation periods to be approx- 
imately equal to the pseudosynchronous period of P ps ~ 1.8 
days derived in §[5] we can constrain i\ and 12'- 



<ii< 11 



and 



2.6° < i 2 < 11 



Wll obtained i orb = 5.52° ± 0.10 by fitting the 
lightcurve's ellipsoidal variation together with radial velocity 
measurements, so the constraints just derived are consistent 
with alignment of rotational and orbital angular momenta, 



Tides act to drive these three inclinations to be mutually 
parallel or antiparallel, so such an alignment once achieved 
is expected to persist indefinitely. In order to simplify the 
analytical formalism as well as reduce the computational 
expense of modeling the observed pulsations, we will adopt 
eq. |TJ as an assumption for the rest of our analysis. 



3 THEORETICAL BACKGROUND 

In this section we review various heterogeneous theoretical 
results that we rely on in later sections. In § |3.1| we sum- 
marize the conventions and definitions used in our analysis. 
In § |3.2| we review the theory of tidally forced adiabatic 
stellar eigenmodes. Later (§ |4), we use this formalism to 
explain qualitative features of the lightcurves of eccentric 
binaries like KOI-54. We also use adiabatic normal modes 
to compute tidal torques (§[5] and Appendix [c| , as well as 
to perform a nonlinear saturation calculation (§ |6.5[ ). How- 
ever, our detailed quantitative modeling of the observations 
of KOI-54 utilizes a nonadiabatic tidally forced stellar oscil- 
lation method that we introduce and employ in § [6] 

In § |3.3| we summarize how perturbed quantities at 
the stellar photosphere, specifically the radial displacement 
and Lagrangian flux perturbation, can be averaged over the 
stellar disk and translated into an observed flux variation. 
Lastly, in § |3.4| we review the traditional approximation, a 
way of simplifying the stellar oscillation equations in the 
presence of rapid rotation. 



3.1 Conventions and definitions 

We label the two stars as per Table [T] consistent with Wll; 
note that the primary/star 1 is taken to be the smaller and 
less massive star. In the following, we focus our analysis 
on star 1, since the results for star 2 are similar. We as- 
sume that both stars' rotational angular momentum vec- 



tors are perpendicular to the orbital plane (§ 2.2 I, and work 
in spherical coordinates (r, 9, <f)) centered on star 1 where 
6 = aligns with the system's orbital angular momentum 
and = points from star 1 to star 2 at periastron. 

We write the stellar separation as D(t) and the true 
anomaly as f(t), so that the position of star 2 is D = 
(D,tt/2, f). We write the semi-major axis as a and the ec- 
centricity as e. The angular position of the observer in these 
coordinates is n = (0 o ,(j) o ), where these angles are related 
to the traditional inclination i and argument of periastron 
to by ( |Arras et aT|2011| ) 

7T 

i = — ~ u mod 2-7T. (2) 



and 



The orbital period [angular frequency] is P or b [f2 or b], 
while a rotation period [angular frequency] is P* [SI*]. The 
effective orbital frequency at periastron is 



S^peri — 



/=0 



1 



1 + e 



(3) 



which is Sip 
ical frequency is 



20. x fi or b for KOI-54. The stellar dynam- 



kJdyn 



GM 



(4) 



^orb 



«2- 



(1) 



which IS CUdyn » 1.1 rad/hr for KOI-54's stars. 
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3.2 Tidal excitation of stellar eigenmodes 

Although we ultimately use an inhomogeneous, nonadia- 
batic code including the Coriolis force to model the pul- 



sations in eccentric binaries (§ 6.2|, the well known normal 



mode formalism provides an excellent qualitative explana- 
tion for many of the features in the lightcurve power spectra 
of systems such as KOI-54. Here we will review the salient 
results of this standard theory; we demonstrate their ap- 
plication to KOI-54 and related systems in §§ |4]-[5] The 
remainder of the paper after § [5] primarily uses our nonadi- 
abatic method described in § |6.2| 

Working exclusively to linear order and operating in the 
coordinates specified in the previous section, we can repre- 
sent the response of star 1 (and similarly for star 2) — all 
oscillation variables such as the radial displacement the 
Lagrangian pressure perturbation Ap, etc. — to a perturbing 
tidal potential by a spatial expansion in normal modes and a 
temporal expansion in orbital harmonics (e.g., |Kumar et ah] 
P95) : 



5X= A n i m k6X nl (r)e- lk ^ bt Y lm {8,(f>). 



(5) 



Here, 2 ^ I < oo and —l^m^l are the spherical harmonic 
quantum numbers, and index the angular expansion; \n\ < 
oo is an eigenfunction's number of radial nodes, and indexes 
the radial expansion^ and |fc| < oo is the orbital Fourier 
harmonic number, which indexes the temporal expansion. 

Each (n,l,m) pair formally corresponds to a distinct 
mode, although the eigenspectrum is degenerate in m since 
for now we are ignoring the influence of rotation on the 
eigenmodes. Each mode has associated with it a set of eigen- 
functions for the various perturbation variables, e.g., £ r , Sp, 
etc., as well as a frequency u„i and a damping rate j n i. For 
stars and modes of interest, 7„; is set by radiative diffusion; 
see the discussion after eq. |l7p . Fig. [T] gives a propagation 
diagram for a stellar model consistent with Wll's mean pa- 
rameters for star 1 (Table [T]). The frequencies of g-modes 
behave asymptotically for n >• and hence w„i <C Wdyn as 
( Christensen-Dalsgaard|2003[ | 

UJ n l ~ LO — , (6) 

n 

where ljo ~ 4 rad/hr for KOI-54's stars. 

The amplitudes A n imk appearing in eq. |5| each rep- 
resent the pairing of a stellar eigenmode with an orbital 
harmonic. Their values are set by the tidal potential, and 
can be expressed analytically: 



T,I 



2 Si CJn I X lm Wlm A„li 

E„i 



(7) 



The coefficients appearing in eq. Q are as follows, 
(i) The tidal parameter e; is given by 



Ah 



Ri 

^Dcr 



(8) 



where 7? p0 ri = a(l— e) is the binary separation at periastron. 

1 Conventionally, n > corresponds to p-modes while n < 
corresponds to g-modes; however, since we are mostly concerned 
with g-modes in this paper, we will report g-mode n values as 
> 0. 
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Figure 1. Propagation diagram for a MESA stellar model consis- 
tent with Wll's mean parameters for star 1 of KOI-54. showing 
the radial profiles of the Brunt- Vaisiila frequency N, the Lamb 
frequency Si = yUJ + l)c 3 /r for I = 2 (where c s is the sound 
speed), and the inverse thermal time 1/ttherm (defined in eq. |16| . 
(AT is dashed where N 2 < 0.) Propagation of g-modes occurs 
where the squared mode frequency ui 2 , is less than both A^ 2 and 



Sf ( Christensen-Dalsgaard 2003|. Nonadiabatic effects become 



important when lo„i ■ ttherm < 27r. Several important frequencies 
for KOI-54 are also plotted, which are defined in § |3.1| The 90th 
harmonic is the largest pulsation observed (Table |2j. 



This factor represents the overall strength of the tide; due 
to its dependence on R\/D pcl \, which is a small number in 
cases of interest, it is often acceptable to consider only / = 2. 

(ii) The linear overlap integral Q n i (Press & Teukolsky 
119771), given by 



1 

MiR\ 
1 

MxR[ 
Ri 



I Kr.nl + (1+ l)£fc,nl) P r dr 



Spni r l+2 dr 



(9) 



GMi 



21 + 1 

47T 



5(j>(Ri 



represents the spatial coupling of the tidal potential to a 
given eigenmode; it is largest for modes with low \n\ and 
hence for eigenfrequencies close to the dynamical frequency 
ojdyn = \J GM\ I R\ , but falls off as a power law for \n\ ^> 0. 
(iii) We define our mode normalization/energy E n i as 



E„i = 2 



J 7llRl 



(£,m + + l)&m) Pr 2 dr, (10) 

LriM! i / Jo 

where £,h is the horizontal displacement ( jChristensen- 
Dalsgaard|2003] . 

(iv) The unit-normalized Hansen coefficients Xf m are the 
Fourier series expansion of the orbital motion (Murray & 
Dermott|1999 1, and are defined implicitly by 

(+1 oo 

e -™/(«) = lL(e)e" W ° bi . (11) 

k — — oo 
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They are related to the traditional Hansen coefficients X u 
by Xf m = Xf m /(l — e) i+1 and satisfy the sum rule 



E 



Xi m {e) 



(12) 



which can be verified u sing eq. (111. (An explicit expression 
for Xi m is given in eq. 



A6 



The Hansen coefficients repre- 
sent the temporal coupling of the tidal potential to a given 



orbital harmonic. They peak near /cp C ak 
fall off exponentially for larger \k\. 
(v) The Lorentzian factor A n i m k is 



2 



mf2p or i/fi or b but 



(13) 



where o^m = fcf2 or b — mS7«, and represents the tempo- 
ral coupling of a given harmonic to a given mode. When 
its corresponding mode/harmonic pair approach resonance, 
i.e. ui n i ~ <Jkm, A n ; m fc can become very large; its maximum, 
for a perfect resonance, is half the simple harmonic oscillator 
quality factor, q n r. A pr = iU rd/2 jni = MjVa/2. 

(vi) Wim is defined in eq. ( A3 1 and represents the angular 



coupling of the tidal potential to a given mode; it is nonzero 
only for mod(Z + m,2) = 0. In particular, W*2,±i = 0, mean- 
ing that I = 2, jmj = 1 modes are not excited by the tidal 
potential. 

(vii) To calculate the quasiadiabatic damping rate 7„; 
within the adiabatic normal mode formalismrl we aver- 
age the product of the thermal diffusivity \ with a mode's 
squared wavenumber k 2 , weighted by the mode energy: 



f rc k\{e + 1(1 + l)e h )pr 2 dr 
fo C (& + l(l + l)e h )pr*dr ' 



where the thermal diffusivity \ i s 

I60-T 3 



3Kp 2 c p 



(14) 



(15) 



The cutoff radius r c is determined by the minimum of the 
mode's outer turning point and the point where w n r£therm = 
2-7T ( Christensen-Dalsgaard 20031, where the thermal time is 

_ pc P T 



gF 



(16) 



When this cutoff is restricted by the mode period intersect- 
ing the thermal time, so that strong nonadiabatic effects 
are present inside the mode's propagation cavity, the mode 
becomes a traveling wave at the surface, and the standing 
wave/adiabatic normal mode approximation becomes less 
realistic. This begins to occur at a frequency (in the rest 
frame of the star) of ~ 50 x fiorb for KOI-54, as can be seen 
in Fig. [l] Fortunately, our calculations involving the normal 
mode formalism (§§ [5] fc |6.5[ ) center primarily on low-order 
modes that are firmly within the standing wave limit. 
The g-mode damping rate scales roughly as 



7ni ~ 7o n 



7o 



up 



(17) 



where we have used the asymptotic g-mode frequency scaling 
from eq. |6]|. In the standing wave regime, i.e. for ui„i < 
50 x fiorb, we find that 70 ~ 1 Myr -1 and s ~ 4. This large 
value for s results from the fact that most of the damping 
occurs at the surface, and the cutoff radius is limited by the 
outer turning point where the mode frequency intersects the 
Lamb frequency. As the mode frequency declines, the cutoff 
radius moves outward toward smaller Lamb frequency and 
stronger damping, as can be seen in Fig. [I] Without this 

behavior of the turning point, we would expect s ~ 2 since 

1 2 2 • 

k oc n m eq. 



(141 



3.3 Observed flux perturbation 

Throughout this work, perturbations to the emitted flux 
AF are understood to be bolometric, i.e. integrated over 
the entire electromagnetic spectrum. We correct for Ke- 
pler's bandpass to first order as follows. We define the band- 
pass correction coefficient j3(T) as the ratio of the bandpass- 
corrected flux perturbation (AF/F)b pc to the bolometric 
perturbation (AF/F), so that 



AF\ 



bpc 



(18) 



We assume Kepler is perfectly sensitive to the wave- 
length band (Ai,A 2 ) = (400,865) nm (|Koch et al.]|2010[), 



and is completely insensitive to all other wavelengths. Then 
P(T) is given to first order by 



P(T) 



J^(dB x /d\nT)d\ 



4AT 



B\d\ 



(19) 



where B\(T) is the Planck function. Using Wll's mean pa- 
rameters for KOI-54 (Table [TJ, we have P(Ti) = 0.81 and 
P(T2) = 0.79. Note that employing /3 alone amounts to ig- 
noring bandpass corrections due to limb darkening. We have 
also ignored the fact that in realistic atmospheres, the per- 
turbed specific intensity depends on perturbations to gravity 
in addition to temperature; this is a small effect, however, 



as shown e.g. in Robinson et al. ( 1982 1 



For completeness, we transcribe several results from 
Pfahl et al. (20081, which allow a radial displacement field 



£ r and a Lagrangian perturbation to the emitted flux AF, 
both evaluated at the stellar surface, to be translated into a 
corresponding disk-averaged observed flux perturbation 8J, 
as seen e.g. by a telescope ( Dziembowski|[l~977 1. While an 
emitted flux perturbation alters the observed flux directly, a 
radial displacement field contributes by perturbing a star's 
cross section|f] 

Given £ r and A_F expanded in spherical harmonics as 



£r = J2 £,r,lm(i)Y lm (e,cj, 
!=0 m=-i 
00 I 

AF = ^ Y, AF lm (t)Y lm (6,< 



(20) 
(21) 



2 We only use this approximate method of calculating damping 
rates when employing the adiabatic normal mode formalism; our 
nonadiabatic method introduced in § |6.2| fully includes radiative 
diffusion. 



3 A horizontal displacement field ^ produces no net effect to 
first order — its influence cancels against perturbations to limb 
darkening, all of which is included in eq. \22\. 
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Table 3. First six disk-integral factors b; and c; from equations 
^23\ and i ]24[ ) for linear Eddington limb darkening, h(p) = 1 + 
(3/2)/.. 



1 ft. 




I 


h 


Q 


1 





3 


1/16 


3/4 


1 17/24 


17/12 


4 


-1/48 


-5/12 


2 13/40 


39/20 


5 


-1/128 


-15/64 



we can translate these into a fractional observed flux varia- 
tion 8J/ J to first order by 



SJ 
J 



= £ £ 



(2& ; - cj) 



R 

AF lm (t) 



F(R) 

where the disk-integral factors are 



bi = pPi(p)h(p) dp 



2 dP 3 ,d 2 P, 
2p — - (p - p 



Y lm {e ,4> \ (22) 
(23) 

h(p) dp, (24) 



dp vr~ f~ / ^2 

Pi in) is a Legendre polynomial, and /i(/Lt) is the limb darken- 
ing function, normalized as J Q ph(p)dp = 1. For simplicity, 
we use Eddington limb darkening for all of our analysis, with 
h(p) — 1 + Sp/2; bi and ci in this case are given in Table [3] 
for sC I < 5. 

Since 5"2,±2 oc sin 2 8 and Yi.o oc (3 cos 2 8 — 1), and since 
KOI-54 has S„ = i?s 5.5° (Table [if, eq. (|22j) shows that 
m = ±2 eigenmodes are a factor of ~ 200 less observable 
than m = modes. It is thus likely that nearly all of the 
observed pulsations in KOI-54 have m = 0; the exceptions 
may be Fl and F2, as we discuss in § |6.6| 

3.4 Rotation in the traditional approximation 

Stellar rotation manifests itself in a star's corotating frame 



as the fictitious centrifugal and Coriolis forces (Unno et al. 
|1989[ ). The centrifugal force directly affects the equilibrium 
structure of a star, which can then consequently affect stel- 
lar oscillations. Its importance, however, is characterized by 
(fi„/u;dyn) 2 , which is ~ 10 -2 for rotation periods and stel- 
lar parameters of interest here (§ [2j. As such we neglect 
rotational modification of the equilibrium stellar structure 
( |Ipser fc Lindblom|1990) . 

The Coriolis force, on the other hand, affects stellar os- 
cillations directly. Given a frequency of oscillation a, the 
influence of the Coriolis force is characterized by the dimen- 
sionless rotation parameter q given by 



2fi*/cx, 



(25) 



where large values of |g| imply that rotation is an important 
effect that must be accounted for. Note that for simplicity 
we assume rigid-body rotation throughout. For the pulsa- 
tions observed in KOI-54's lightcurve (Table assuming 
both stars rotate at near the pseudosynchronous rotation 
period P ps ~ 1.8 days discussed in § |5.1| and that m = 
(justified in the previous section), q ranges from 0.5 for 
k — cr/Sl or b = 90 to 1.5 for k — 30. Thus lower harmonics 



fall in the nonperturbative rotation regime, where rotation 
is a critical effect that must be fully included. 



The "traditional approximation" ( |Chapman fc Lindzen 



1970 1 greatly simplifies the required analysis when the Cori- 
olis force is included in the momentum equation. In the case 
of g-modes, it is applicable for 



1 > 



R 



\N\ 



(26) 



where H p — pg/p is the pressure scale height; outside of 
the convective cores of models we are concerned with in this 



work (where g-modes are evanescent anyway), eq. ( 26 1 is well 
satisfied whenever rotation is significant. Here we will give 
a brief overview of the traditional approximation; we refer 
to Bildsten et al. ( 1996 1 for a more thorough discussion. 

The traditional approximation changes the angular 
Laplacian, which occurs when deriving the nonrotating stel- 
lar oscillation equations, into the Laplace tidal operator 
L^. (Without the traditional approximation, the oscillation 
equations for a rotating star are generally not separable.) It 
is thus necessary to perform the polar expansion of oscilla- 
tion variables in eigenfunctions of L^, known as the Hough 
functions H^ m (p) (where p = cos 6), rather than associated 
Legendre functions; the azimuthal expansion is still in e !m< *. 
The eigenvalues of are denoted A, and depend on m, the 
azimuthal wavenumber. In the limit that q — > 0, the Hough 
functions become ordinary (appropriately normalized) asso- 
ciated Legendre functions, while A — >■ 1(1 + 1). 

We present the inhomogeneous, tidally driven stellar os- 
cillation equations in the traditional approximation in Ap- 
pendix |A"2) The principal difference relative to the standard 
stellar oscillation equations is that terms involving 1(1 + 1) ei- 
ther are approximated to zero, or have 1(1 + 1) — > A. This re- 
placement changes the effective angular wavenumber. E.g., 
since the primary A for m — increases with increasing ro- 
tation, fast rotation leads to increased damping of m — 
g-modes at fixed frequency, as discussed in § |6.3| 

For strong rotation, \q\ > 1, the Hough eigenvalues 
A can be both positive and negative. The case of A > 
produces rotationally modified traditional g-modes, which 
evanesce for cos 2 8 > 1/g 2 . (Rossby waves or r-modes are 
also confined near the equator and have a small positive 
value of A.) Instead, for A < 0, polar modes are produced 
that propagate near the poles for cos 2 8 > 1/q 2 , but evanesce 
radially from the surface since they have an imaginary Lamb 
frequency S\ = A 1 ^ 2 c s /r (as explained further in Fig. [lj. 



4 QUALITATIVE DISCUSSION OF TIDAL 
ASTEROSEISMOLOGY 

It is helpful conceptually to divide the tidal response of a star 
into two components, the equilibrium tide and the dynamical 
tide ( Zahn|1975[ ) . Note that in this section we will again use 
the normal mode formalism described in § |3.2| even though 
our subsequent more detailed modeling of KOI-54 uses the 
inhomogeneous, nonadiabatic formalism introduced in §|6.2| 



4.1 Equilibrium tide 

The equilibrium tide is the "static" response of a star to 
a perturbing tidal potential, i.e., the large-scale prolation 
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due to differential gravity from a companion. In terms of 
lightcurves, the equilibrium tide corresponds to ellipsoidal 
variability (along with the irradiation component of this ef- 
fect discussed in Appendix Bl I. In the case of an eccentric 



binary this manifests itself as a large variation in the ob- 
served flux from the binary during periastron. KOI-54's 
equilibrium tide was successfully modeled in Wll, enabling 
precise constraints to be placed on various stellar and orbital 
parameters (Table [T]). 

In terms of the normal mode formalism developed in 
§ |3,2| the equilibrium tide corresponds to the amplitudes 
from eq. |7f tied to large overlaps Q„i and large Hansen 
coefficients Xf m \ in other words, to pairings of low-jnj modes 
with low-|fc| orbital harmonics. The Lorentzian factor A„i m k 
is typically ~ 1 for the equilibrium tide since it is not a 
resonant phenomenon. 

In practice, however, it is much simpler and more con- 
venient to use other mathematical formalisms to model the 
equilibrium tide, like taking the zero-frequency stellar re- 
sponse as in Appendix |B2| or filling Roche potentials as in 
Wll's simulations. We show in § |6.1| that our simple ana- 
lytical treatment of the equilibrium tide verifies the results 
from the sophisticated simulation code employed in Wll. 



4.2 Dynamical tide 

The dynamical tide, on the other hand, corresponds to reso- 
nantly excited pulsations with frequencies equal to harmon- 
ics of the orbital frequency, kQ OI h- Wll observed at least 
21 such harmonics (Table Ejl. 

For a circular orbit, the tidal potential has all its power 
in the k = ±2 orbital harmonics; in this case the only 
modes that can be resonantly excited are those with fre- 
quencies close to twice the Doppler-shifted orbital frequency: 
oj n i w 2|fi or b — 0*|; this is typically only a single mode. 
This corresponds to the fact that the Hansen coefficients 
from eq. ( |11[ ) become a Kronecker delta at zero eccentricity: 
^im(O) = &m- However, for a highly eccentric orbit, the 
distribution of power in the Hansen coefficients, and hence 
the stellar response, can be much broader; as a result a wide 
array of different harmonics can be excited, allowing for a 
rich pulsation spectrum. 

Mode excitation due to a tidal harmonic fcf2 or b is modu- 
lated by the Doppler-shifted frequency a^ m = fcf2 or b — mQ„. 
However, the frequencies at which modes are observed to 
oscillate, viewed from an inertial frame, are indeed pure 
harmonics of the orbital frequency, fcf2 or b|^] We demon- 
strate this mathematically in Appendix [A] intuitively, al- 
though a driving frequency experiences a Doppler shift upon 
switching to a star's corotating frame, the star's response is 
then Doppler shifted back upon observation from an inertial 
frame. In general, any time a linear system is driven at a 
particular frequency, it then also oscillates at that frequency, 
with its internal structure reflected only in the oscillation's 
amplitude and phase. 

Whether a given mode is excited to a large amplitude 
is contingent on several conditions — essentially all the terms 



in eq. Q. First, the overall strength of the tide, and hence 
the magnitude of observed flux variations, is determined by 
the tidal factor ei from eq. Q. The dominant multipole 
order is I = 2, so we have £2 = {Mi/Mi){Rx/ 'D peT if where 
D pcl i = a(l — e) is the binary separation at periastron, and 
we are focusing our analysis on star 1. For KOI-54, £2 — 
4 x 1CT 3 for both stars. 

Next, the strength of a mode's resonant temporal cou- 
pling to the tidal potential is given by the Lorentzian factor 
A n ; m fe in eq. (13 1. Since this factor is set by how close a 



mode's frequency is to the nearest orbital harmonic, its ef- 
fect is intrinsically random. The degree of resonance has an 
enormous effect on a mode's contribution to the observed 
flux perturbation, meaning that modeling the dynamical 
tide amounts on some level to adjusting stellar and system 
parameters in order to align eigenfrequencies against orbital 
harmonics so that the array of Lorentzian factors conspire 
to reproduce observational data. 

Moreover, given a single observed pulsation amplitude 
together with theoretical knowledge of the likely responsible 
mode, i.e. the first four factors in eq. Q, equating theoret- 
ical and observed pulsation amplitudes in principle yields 
direct determination of the mode's eigenfrequency, indepen- 
dently of the degree of resonance. This line of reasoning 
of course neglects the considerable theoretical uncertainties 
present, but serves to illustrate tidal asteroseismology's po- 
tential to constrain stellar parameters. 

Despite the inherent unpredictability, a lightcurve's 
Fourier spectrum is still subject to restrictions imposed pri- 
marily by the remaining two factors in eq. These terms, 
the linear overlap integral Q n i and the unit-norma lized 
Hansen coefficient Xi m {e) (respectively equations [9] and 11 1, 
restrict the range in k over which harmonics can be excited; 
Fig. [2] shows profiles of both. As discussed in § |3.2[ Q n i 
peaks for modes with frequencies near the dynamical fre- 
quency of the star ojd yn and falls off as a power law in fre- 
quency, whereas Xf m peaks for harmonics near mQ. pe ri/0, OI b 
and falls off for higher k: 

Qni CC UJ n l < ^dyn (27) 

X /m (e) oc exp(-k/r) \k\ > mtt pcr i /Q orh . (28) 
The power-law index p is 11/6 for g-modes in stars with a 



convective core and a radiative envelope or vice versa ( Zahn 



19701, and for KOI-54's eccentricity and I = 2 we find r 



15. 

As a result, modes that can be excited are those with 
frequencies in the intervening region between the peaks of 
Q n i and Xf m , i.e., 



|m|Jl P eri < U n i < Wdyn- 



(29) 



4 |Welsh et al.| l |201l| | incorrectly attributed nonharmonic pulsa- 
tions to rotational splitting; we return to nonharmonic pulsations 
in S|675l 



This is a necessary but not sufficient condition; Fig.[3]shows 
the product QniX} s m (e) at various eccentricities with stel- 
lar parameters as well as the periastron distance D pcr j fixed 
to the mean values in Wll, and hence with fixed tidal pa- 
rameter ei (but consequently allowing the orbital period to 
vary). Although a chance close resonance can yield a large 
Lorentzian factor A n i m fc, excitation of modes far from the 
peak of QnlXim becomes less and less likely, since this quan- 
tity falls off sharply, especially towards larger \k\. 

There are two other constraints on the range of har- 
monic pulsations that can be excited. First, the eigenmodc 
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Figure 2. The linear overlap integral Q n i and unit-normalized 
Hansen coefficients X* (e) as a function of k = ui/il or -^ for 1 = 2 
and m = 0,2 (note the identity Xf = X^^ m ). Stellar and 
orbital parameters are fixed to Wll's mean values for star 1; 
in particular, e = 0.8342. Each point on the curve for Q repre- 
sents a normal mode with frequency u) n i ~ fcf2 or t,, thus neglecting 
any Doppler shift due to rotation. (See § |6.3| for a discussion of 
the influence of Doppler shifts.) Modes of a given m can be ex- 
cited near where the overlap curve intersects the Hansen curve, 
at log 10 (fc) m 2.0 in this plot. 



density for g-modes scales asymptotically as 

Wdyn 



dn 
dk 



J_ 



(30) 



which shows that fewer modes exist at higher k. This can be 
seen by the spacing of points (which denote normal modes) 
in Figures [2] and [3] as well as by the spacing of peaks in 
Fig. [6] This further limits the number of harmonics that 
can be excited at large k, in addition to the exponential 
decay of the Hansen coefficients discussed earlier, and thus 
effectively shifts the curves in Fig. [3] toward lower k. 

In addition, the Lorentzian factor A n i m k is attenuated 
by mode damping y n i, which is set by radiative diffusion 
for high-order g-modes. Damping becomes larger with de- 
creasing g-mode frequency due to increasing wavenumber; 
an asymptotic scaling is given in eq. (17 1. Because the 



Lorentzian response is proportional to j~ [ at perfect reso- 
nance, the amplitudes of lower-frequency modes/harmonics 
are diminished by increased damping, in addition to the 
power-law decay of the tidal overlap. This effect is criti- 
cal for understanding the influence of rotation on lightcurve 
power spectra, as we investigate in § |6.3| 



4.3 Pulsation phases 

Pulsation phases in eccentric binaries are essential informa- 
tion which should be fully modeled, in addition to the pul- 
sation amplitudes reported in Wll. For simplicity, we focus 
on a particular harmonic amplitude A„imk from equations 
|5| and Q and assume it results from a close resonance so 
that w (jfc m = fcfiorb — mf2*, assuming without loss of 
generality that a km > 0. We can then evaluate its phase 
ipnimk relative to periastron, modulo 7r (since we are tem- 
porarily ignoring the real part of the amplitude, which could 



lo §io (<2„2 X 2o( e )) 




Figure 3. Plot of Q n iXl L(e) as a function of k = ui/Q^^ for 
several eccentricities with I = 2, m = 0, and all stellar parameters 
as well as the periastron separation D peIl fixed to Wll's mean 
values for star 1 (Table [TJ. Fixing D per i fixes the tidal factor £; 
from eq. (JsJ and hence the overall strength of the tide (although 
the orbital period consequently varies). Each point represents 
a normal mode; Hansen coefficients are evaluated with k given 
by the integer nearest to aj n ; /f2 or t, for each eigenmode. Solid 
vertical lines denote the n = 7 g-mode (the g-mode with 7 radial 
nodes), while dashed vertical lines are n = 14. Note that the finer 
cigcnfrequency spacing at small k allows for larger amplitudes, 
which is not accounted for in this plot; including this effect would 
shift the curves toward lower k. 



introduce a minus sign), 

Ipnimk = arg( J 4„; m fc) 

= 7r/2 — arctan 



r, 



\ l^nlOkm 

7r/2 — arctan(<5o;/7 n ;) mod 7T, 



mod 7T (31) 



where 8ui = uj n i — (Tkm is the detuning frequency. 

For a near-perfect resonance, where 5uj < 7„;, ipnimk 
approaches it/ 2 (modulo it). However, if eigenmode damp- 
ing rates are much smaller than the orbital frequency, then 
this intrinsic phase should instead be near 0. This is the case 
for KOI-54, where £l rb/7ni > 10 3 for modes of interest. In- 
deed, theoretically modeling the largest-amplitude 90th and 
91st harmonics of KOI-54 assuming they are m = modes 
requires only |5w|/7ni ~ 20, so that even these phases should 
be within ~ 1% of zero (modulo n). 

The phase of the corresponding observed harmonic flux 
perturbation can be obtained from eq. (31 1 by further in- 



cluding the phase of the spherical harmonic factor in the 
disk-averaging formula, eq. ( |22[ ): 

(32) 



arg((SJ fc /J) = i/)„imk + m4> B 



Summing over the complex conjugate pair, the observed 
time dependence is then + m4>o)], where 

t = corresponds to periastron. However, the sign of k in 
this formalism is unknown; equivalently, whether the pul- 



sation is prograde or retrograde (Appendix CI I cannot be 
determined in this way. Thus if the observed pulsation's 
(cosine) phase is S, the comparison to make is 



8 = ±{i>„i m k + m<f> ) mod n. 



(33) 
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Nonetheless, since we have argued that i/'nimfc ~ 0, this be- 
comes 



as derived in Hut 



±rruf> mod -n. 



(34) 



Given determination of (f> (related to the argument of 
periastron uj by eq. [2J by modeling of RV data or ellipsoidal 
variation, the phase of a resonant harmonic thus directly 
gives the mode's value of \m\ (which is very likely or 2 for 
tidally excited modes, since I = 2 dominates). For KOI-54, 
phase information on harmonics 90 and 91 would thus de- 
termine whether they result from resonance locks, as dis- 
cussed in the next section, or are simply chance resonances. 
Furthermore, knowing \m\ allows m<f) to be removed from 
eq. (33 1, yielding the pulsation's damping-to-detuning ratio. 

However, the preceding treatment is only valid if the 
eigenfunction itself has a small phase: although eigenfunc- 
tions are purely real for adiabatic normal modes, local 
phases are introduced in a fully nonadiabatic calculation, 



6.2 Thus equations (32 I - (34 1 are only applicable 



in the standing wave regime, where the imaginary part of 
the flux perturbation is small relative to the real part. In 
the traveling wave regime, the local wave phase near the 
surface becomes significant, and can overwhelm the contri- 
bution from global damping; see § |6.2| For KOF54, this 
corresponds to |fc| below ~ 30, although this depends on the 
rotation rate (8 16.31). 



5 ROTATIONAL SYNCHRONIZATION IN 
KOI-54 

Here we will discuss a priori theoretical expectations for 
KOI-54's stars' rotation. Later, in § |6.3| we will compare 
the results derived here with constraints imposed by the 
observed pulsation spectrum. 



5.1 Pseudosynchronization 

In binary systems, the influence of tides causes each compo- 
nent of the binary to eventually synchronize its rotational 
and orbital motions, just as with Earth's moon. Tides 
also circularize orbits, sending e — > 0, but the circulariza- 
tion timescale t c i IC is much greater than the synchronization 
timescale i sync ; their ratio is roughly given by the ratio of 



orbital to rotational angular momenta: 



Lorb 

(if 



(35) 



(1 



where 7* is the stellar moment of inertia, fi is the reduced 
mass, and we have assumed for simplicity that the stars 



rotate at the periastron frequency r2 pcr i (§ 3.1 ). For KOI-54, 
this ratio is ~ 10 3 . 

Due to the disparity of these timescales, a star in an ec- 
centric binary will first synchronize to a pseudo synchronous 
period P ps , defined as a rotation period such that no aver- 
age tidal torque is exerted on either star over a sufficiently 
long timescale. If only the torque due to the equilibrium 
tide is used, and thus eigenmode resonances are neglected, 
then only one unique pseudosynchronous period exists, Pp S r , 



for KOI-54 is (eq. 



(19811 and employed in Wll. Its value 



C13I 



(2.53 ± 0.01) days. 



Inclusion of eigenmode resonances, however, compli- 
cates the situation. Fig. [4] shows the secular tidal torque 
(averaged over one rotation period) for star 1 of KOI-54 
plotted as a function of rotation frequency/period includ- 
ing contributions from both the equilibrium and dynamical 
tide. Although the general torque profile tends to zero at 
Pps , numerous other roots exist (displayed as vertical lines), 
where the torque due to a single resonantly excited eigen- 
mode of the dynamical tide cancels against that due to the 
equilibrium tide. To produce this plot, we directly evaluated 



the secular tidal torque (Appendix CI I using an expansion 
over the quadrupolar adiabatic normal modes of a MESA 



stellar model (Paxton et al. 20111 with parameters set by 



Wll's mean values for star 1 of KOI-54 (Table [T}. In our 
calculation we include both radiative (§ 3.2| and turbulent 
convective damping ( Willems et al.|2010 |, but neglect rota- 
tional modification of the eigenmodes. 

Next, of the many zeroes of the secular torque avail- 
able, which are applicable? Continuing with the assump- 
tion that KOI-54's stars were born with rotation periods 
of Pbirth ~ 1 day (§ 2.1 1, with the same orientation as the 
orbital motion, one might naively posit that the first zero 
encountered by each star should constitute a pseudosyn- 
chronous period — it is an ostensibly stable spin state since 
small changes to either the stellar eigenmodes (via stellar 
evolution) or the orbital parameters (via circularization and 
orbital decay) induce a restoring torque. This is the basic 
idea behind a resonance lock ( |Witte fc Savonije|1999| |. 

However, this conception of resonance locking neglects 
two important factors. First, although the dynamical and 
equilibrium tidal torques may cancel, their energy deposition 
rates do not (in general); see Appendix CI Thus during a 



resonance lock the orbital frequency must continue to evolve, 
allowing other modes to come into resonance, potentially ca- 
pable of breaking the lock. Second, as shown by |Fuller fc Lai| 
( |2011[ ), it is necessary that the orbital frequency not evolve 
so quickly that the restoring torque mentioned earlier be in- 
sufficient to maintain the resonance lock. This restricts the 
range of modes capable of resonantly locking, introducing an 
upper bound on their inertial-frame frequencies and hence 
their orbital harmonic numbers (values of k in our notation). 

Consequently, pseudosynchronization is in reality a 
complicated and dynamical process, consisting of a chain 
of resonance locks persisting until eventually e — > and 
P» = P or b. Such resonance lock chains were studied in 



much greater detail by Witte & Savonije (1999) for eccen- 



tric binaries broadly similar to KOI-54. As a result of the 
inherent complexity, a full simulation of KOI-54's orbital 
and rotational evolution is required in order to address the 
phenomenon of resonance locking and to derive theoretical 
predictions for the stars' spins. To perform such simulations, 
we again expanded the secular tidal torque and energy de- 



position rate over normal modes (detailed in Appendix CI I 



using two MESA stellar models consistent with Wll's mean 
parameters for KOI-54's two stars. We then numerically in- 



tegrated the orbital evolution equations (Witte & Savonije 
|1999[ ) assuming rigid-body rotation. We did not include the 
Coriolis force, nor did we address whether the eigenmode 
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Figure 4. Plot of the secular tidal synchronization torque r as 
a function of the rotational frequency /„ = 1/P* for star 1, us- 
ing Wll's parameters (Table Top panel: Black indicates 
a positive torque (meaning increasing stellar spin), while red in- 
dicates a negative torque (decreasing stellar spin). The over- 
all profile of r goes to zero at Hut's value for the nonresonant 
pseudosynchronous period, Pp S r = 2.53 days for KOI-54, but 
eigenmode resonances create many additional zeroes, displayed 
as light gray vertical lines. (Determination of zeroes in this plot 
is limited by its grid resolution. Many more exist; see bottom 
panel.) Bottom panel: Zoom in, showing how narrow resonance 
spikes can cause the otherwise negative torque to become zero. 
The torque pattern roughly repeats at half the orbital frequency, 
/orb/2 = 0.012 day -1 , due to the fact that eigenmodes resonate 
with Doppler-shifted driving frequencies <r fem = kQ, Q ^ — mil,, 
with m = ±2 and k S Z. 



appears independently of the initial rotation rates used; in 
other words, it is an attractor. 

As described above, when a star is locked in resonance, 
it is the torque from a single highly resonant eigenmode 
that acts to oppose the equilibrium tide's nonzero torque. 
Such a high-amplitude mode should be easily observable. 
At first glance, this line of reasoning seems to provide a 
natural explanation for the presence of the large-amplitude 
90th and 91st observed harmonics in KOI-54 (Fl and F2 
from Table[5|, namely that each is the photometric signature 
of the highly resonant eigenmode that produces a resonance 
lock for its respective star. There are several problems with 
this idea, however, which we elucidate in 8 16.61 



5.2 Synchronization timescale 

Where between the stars' putative birth rotation periods, 
fbirth = 1.0 day, and the pseudosynchronous period from 
our simulations, P ps ~ 1.8 days, do we a priori expect the 
rotation periods of KOI-54's stars to fall? To this end, we 
can roughly estimate the synchronization timescale t Bync by 
integrating IQ, — t(S1») to find 



ts 



d f2„ 



(36) 



and 



where / is the stellar moment of inertia, Q ps = 2ir/P ps , 
the tidal torque r can as before be calculated as a function 
of the spin frequency il* using an expansion over normal 



modes (Appendix CI I 



Using this approximation, we find t sync ~ 80 Myr, 
which is less than the inferred system age of t agc ~ 200 Myr. 
This is consistent with our orbital evolution simulations 
(§ |5.1[ ). Although stellar evolution was ignored in this calcu- 
lation, a rough estimate of its effect can be made using only 
the fact that t Bync scales as R~ 3 (since the torque scales as 
R 5 while the moment of inertia scales as R 2 ). Given that 
both stars had 10% smaller radii at ZAMS (indicated by our 
modeling), this would lead to only at most a 3 x 10% = 30% 
increase in t aync . Furthermore, both stars had much larger 
radii before reaching the main sequence, which would imply 
an even shorter synchronization time. Lastly, an important 
effect that arises when rotation is fully included is the exis- 
tence of retrograde r-modes, which would also enhance the 
rate of stellar spindown ( Witte fc Savonije|19 99). Thus the 
inequality 



> t B 



seems to be well satisfied, and we expect that both stars' 
rotation periods should be close to the value of P ps ~ 1.8 
days from § |5.1| 



amplitudes required to produce the various resonance locks 



that arise are stable to nonlinear processes (§ 6.5 I 



Our simulations indicate that both stars should have 
reached pseudosynchronization states with rotation periods 



in more detail in § |5.2| These periods are ~ 30% faster than 
Hut's value of Pp S r = 2.5 days. The pseudosynchronization 
mechanism that operates is stochastic in nature, in which 
the dynamical tide's prograde resonance locks balance the 
equilibrium tide in a temporally averaged sense. This result 



6 RESULTS 

6.1 Ellipsoidal variation 

Fig. |5]a shows our simple model of KOI-54's ellipsoidal 
variation; we adopted the best-fit parameters from Wll's 
modeling (Table [T]) to produce our lightcurve. Our irra- 
diation (Appendix |B1| blue dashed line) and equilibrium 
tide (Appendix |B2| red dotted line) models are larger than 
Wll's results by 24% and 14% respectively. The shapes of 
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Figure 5. Simple analytic model of ellipsoidal variability detailed in Appendix[B] including both the equilibrium tide (red dotted lines) 
and "reflection" /irradiation (blue dashed lines) components of the lightcurve. We used the best-fit parameters from Wll in all three 
panels (Table [TJ, except that we show two examples of edge-on orientations (ignoring the possibility of eclipses) in (b) & (c), effectively 
presenting KOI-54's lightcurve as it would be observed from different angles. (We used u> = (80°, 20°) for (b, c) in order to demonstrate 
the asymmetric lightcurves possible depending on the binary's orientation.) Panel (a) reproduces Wll's modeling and the data for 
KOI-54 to ~ 20%. Our analytic model is easily applicable to many other systems. In § |6.3| we show that the dynamical tidal response, 
ignored here, may be larger than that due to the equilibrium tide for edge-on systems. 



both curves are, however, essentially indistinguishable from 
Wll's much more detailed calculations. 

We attribute the small difference between our results 
and those of Wll to our simple model of the bandpass cor- 
rection (eq. 19 1 which ignores bandpass variations due to 
limb darkening. Such details could easily be incorporated 
into our analytical formalism, however, by introducing a 
wavelength-dependent limb darkening function h\ (/i) in the 
disk integrals in equations ( |23[ ) and ( |24| ) (Ro binson et al.] 
|1982[ ). We thus believe that the models provided in Ap- 
pendix [B] should be quite useful for modeling other systems 
like KOI-54, due in particular to their analytic simplicity. 

We also show in Fig.[5]b & c what KOI-54's equilibrium 
tide and irradiation would look like for two edge-on orien- 
tations, demonstrating the more complicated, asymmetric 
lightcurve morphologies possible in eccentric binaries (see 
also the earlier work by Kumar et al.|1995 l. Future searches 
for eccentric binaries using Kepler and other telescopes with 
high- precision photometry should allow for the wide range 
of lightcurve shapes shown in Fig. [5] We note, however, 
that that the dynamical tidal response, ignored in this sec- 
tion, may be larger than that due to the equilibrium tide for 
edge-on systems, as we show in § |6.3| 



6.2 Nonadiabatic inhomogeneous method 

Thus far our theoretical results have primarily utilized 
the tidally forced adiabatic normal mode formalism. Al- 
though this framework provides excellent intuition for the 
key physics in eccentric binaries, it is insufficient for produc- 
ing detailed theoretical lightcurves, since this necessitates 
tracking a star's tidal response all the way to the photo- 
sphere where nonadiabatic effects are critical. To account for 
this, we employ the nonadiabatic inhomogeneous formalism 
originally used by Pfahl et al. ( 2008[ ) (Appendix |A1[ ), which 
we have extended to account for rotation in the traditional 



Rather than decompose the response of the star into 
normal modes, the inhomogeneous method directly solves 
for the full linear response of the star to an external tidal 
force produced by a companion at a given forcing frequency. 
Given a stellar model, an orbital period, a set of orbital 
harmonics to act as driving frequencies, and a rigid-body 
rotation period, we solve the numerical problem described 
in Appendix |A2| for each star. This determines the various 
physical perturbation variables of the star as a function of 
radius, such as the radial displacement and the flux pertur- 
bation. For stars of interest we can safely ignore perturba- 
tions to the convective flux, so the only nonadiabatic effect 
is that produced by radiative diffusion. 

Fig. [6] shows the surface radial displacement £ r and La- 
grangian emitted flux perturbation A_F computed on a fine 
frequency grid, temporarily ignoring rotation; normal mode 
frequencies correspond to the resonant peaks in these curves. 
The surface radial displacement should approach its equi- 
librium tide value as the driving frequency tends to zero. 
Quantitatively, we find that this is true for orbital harmonics 
k < 30; note that in the units employed in Fig. [6j this equi- 
librium tide value for £ r is (£r/-R)phot = 1 (Appendix B2|. 



approximation (Appendix A2 1 



The surface flux perturbation shown in Fig. [6] on the 
other hand, more clearly demonstrates the three qualita- 
tively different regimes possible at the surface. First, the 
weakly damped standing wave regime, k > 30, is charac- 
terized by strong eigenmode resonances and all perturba- 
tion variables having small imaginary parts. In Fig. [TJ this 
corresponds to the outer turning point, where the mode fre- 
quency intersects the Lamb frequency, lying inside the point 
where the mode frequency becomes comparable to the ther- 
mal frequency, so that the mode becomes evanescent before 
it becomes strongly nonadiabatic. 

Next, the traveling wave regime, 5 < k < 30, arises 
when modes instead propagate beyond where the mode and 
thermal frequencies become comparable, leading to rapid 
radiative diffusion near the surface. In the traveling wave 
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Figure 6. Amplitude of both the radial displacement £ r and Lagrangian flux perturbation A_F evaluated at the photosphere as a function 
of k = aj/Q OT b for a MESA model of star f , using Wlf 's best-fit parameters (TableJTJ. As a result of its larger amplitude, the Lagrangian 
flux perturbation has a much larger effect than the surface displacement on observed flux variations. Rotation is not included in this 
calculation. The radial dependence of the tidal potential is set to U = —(GMi/Ri)(r/Ri) 1 with I = 2, which determines the vertical 
scale. Normal mode eigenfrequencies correspond to peaks in the curves. See the text for a discussion of the different regimes present in 
the stellar response for different forcing frequencies. 



limit, resonances become severely attenuated as waves are 
increasingly unable to reflect at the surface, and all pertur- 
bation variables have comparable real and imaginary parts 
(not including their equilibrium tide values). 

Lastly, just as with the radial displacement, the flux 
perturbation also asymptotes to its overdamped equilib- 
rium tide/von Zeipel value of (AF/F) c p ^ t A = -(I + 
2)(£ r /-R)p') ! lot (Appendix B2| in the low-frequency limit, 
which is lAF/Fj^f = | - I - 2| = 4 in Fig. [els units. 
Quantitatively, however, this only occurs for k < fx At first 
glance, this suggests that the equilibrium tide modeling of 
KOI-54 in Wll and Fig. [5] is invalid, since the equilibrium 
tide in KOI-54 has orbital power out to at least k ~ 30 (as 
can be seen e.g. in the plot of the Hansen coefficients for 
KOI-54's eccentricity in Fig. [2). 

Fortunately, as we describe in the next section, includ- 
ing rotation with a face-on inclination effectively stretches 
the graph in Fig. [6] towards higher k. E.g., for P* = 2.0 
days, we find the equilibrium tide/von Zeipel approxima- 
tion to hold for k < 30, justifying the simplifications used 
in Wll and Appendix |B2[ although this may not apply for 
edge-on systems. 



6.3 Effect of rotation on the dynamical tidal 
response 

The most important effects of rotation in the context of 
tidal asteroseismology can be seen in Fig. [7] Here we show 
the predicted flux perturbation for KOI-54 as a function of 
orbital harmonic k for four different rotation periods, hav- 
ing subtracted the equilibrium tide (Appendix |B2[ ) to focus 
on resonant effects. [^] In Fig. [TJa we use KOI-54's face-on 
inclination of i = 5.5°, while in Fig. [7]b we use an incli- 
nation of i = 90° to illustrate how a system like KOI-54 
would appear if seen edge on; all other parameters are fixed 
to those from Wll's modeling (and are thus not intended 
to quantitatively reproduce the data; see Fig. [8] for an op- 
timized model). The details of which specific higher har- 
monics have the most power vary as rotation changes mode 
eigenfrequencies, moving eigenmodes into and out of reso- 



5 We assume that rotation is in the same sense as orbital motion 
throughout this section. 



nance. Nonetheless, several qualitative features can be ob- 
served. 

For KOT54's actual face-on orientation, as in Fig. [7]a, 
rotation tends to suppress power in lower harmonics. This 
can be understood as follows. Primarily m = modes are 
observable face on (§ |3.3[ ). At fixed driving frequency a, 
as the stellar rotation frequency f2«, and hence the Coriolis 
parameter q — 2f2,/cr, increases in magnitude, m — g- 
modes become progressively confined to the stellar equator 



(§ 3.4 1. As a result, these rotationally modified modes angu- 
larly couple more weakly to the tidal potential, diminishing 
their intrinsic amplitudes. Moreover, equatorial compres- 
sion also corresponds to an increase in the effective multi- 
pole I, where I ~ y/X, and A is a Hough eigenvalue from 
§[373] (e.g., Fig. 2 of |Bildsten et al.|p9%]) Consequently, 
since g-modes asymptotically satisfy eq. (|6j), the number of 
radial nodes n must increase commensurately. Larger n in- 
creases the radial wavenumber, which enhances the damping 
rate, further suppressing the resonant response of the modes 
and hence their contribution to the observed flux variation. 
This effectively corresponds to extending the highly damped 
traveling wave regime toward higher k in Fig. [6] 

As described in § |3.4[ when the magnitude of the Cori- 
olis parameter becomes greater than unity, a new branch 
of eigenmodes develops with negative Hough eigenvalues, 
A < 0. These modes are confined to the stellar poles rather 
than the equator (Lindzcn 1966). They also have an imag- 
inary Lamb frequency, so that they are radially evanescent 
(explained further in Fig.[TJ), and couple weakly to the tidal 
potential. We found negative-A modes to produce only a 
small contribution to the stellar response, which increased 
with increasing rotation rate but which was roughly con- 
stant as a function of forcing frequency, thus mimicking the 
equilibrium tide. The role of these modes in the context 
of tidal asteroseismology should be investigated further, but 
for now we have neglected them in Fig. [7] 

For edge-on orbits, as in Fig. [7|b, the situation is more 
complicated, and there are high-amplitude pulsations ob- 
servable at all rotation periods. First, m — modes very 
weakly affect edge-on lightcurves, since their Hansen coef- 
ficients (which peak at k = 0) do not intersect with the 
linear overlap integrals as strongly as for m — 2 modes (ex- 
plained further in § |4.2| and shown in Fig. 121. Similarly, 
modes with m — — 2 have Hansen coefficients which peak 
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Figure 7. Influence of rotation on the lightcurve temporal power spectra of eccentric binaries. The eight leftmost plots show theoretical 
power spectra using fiducial stellar models and orbital parameters consistent with Wll's mean values (Table [TJ, each as a function of 
orbital harmonic k. The top row uses KOI-54's known face-on inclination of i = 5.5°; the bottom row instead uses an edge-on inclination 
of i = 90° , showing how a system like KOI-54 would appear if viewed edge on. The vertical axis has different scales in the two rows. The 
effects of rotation on the stellar pulsations are included using the traditional approximation (§ |3.4[ |. The equilibrium tide (Appendix |B2[ l 
has been subtracted in order to focus on resonant effects. We have not included negative-A modes, as discussed in the text. (Note that 
parameters here have not been optimized to reproduce KOI-54's lightcurve; see Fig.[8]for such a model.) The four leftmost columns show 
four different rigid-body rotation periods, (a) For a perfectly face-on orientation, only m = modes can be observed (§ |3.3| ; however, 
larger rotation rates lead to equatorial compression of lower-frequency m = g-modes, which increases the effective I and hence enhances 
the dissipation, leading to attenuated amplitudes. The rightmost panel shows the data for KOI-54, which is qualitatively most consistent 
with shorter rotation periods of 1.0 - 2.0 days, comparable to the pseudosynchronous period of ~ 1.8 days calculated in § [5] (b) For 
edge-on systems, a similar argument regarding rotational suppression of mode amplitudes applies, but instead near the Doppler-shifted 
harmonic k = 2P or t,/P«, rather than for k near zero for the m = modes observable face on (see text for details). Comparison of the four 
left panels to the rightmost panel, which shows our simple analytical equilibrium tide model's harmonic decomposition (Appendix |B2| , 
demonstrates that the dynamical tide can dominate the lightcurve in edge-on systems. 



near — 2f2p cr ; /f2 or b and sxg very small for k ^ oQ Thus re- 
gardless of rotation, only m = +2 modes make significant 
lightcurve contributions. 

Within the m = +2 modes, there are two regimes 
to consider: prograde modes excited by harmonics k > 
2f2„/f2 or b and retrograde modes with k < 2f2,/f2 or b (see 



Appendix CI I. Prograde modes at a given corotating frame 
frequency are Doppler shifted toward large k, whereas the 
Hansen coefficients peak near 2f2 pcr i/f2 or b, so their contri- 
bution to lightcurves is marginalized for fast rotation. 

Retrograde, m = +2 g-modes with small corotating- 
frame frequencies atm = fcf2 or b — mf2* are subject to the 
same effect described earlier in the face-on case: they are 
suppressed by fast rotation due to weaker angular tidal cou- 
pling and stronger damping. The difference, however, is that 
although small driving frequencies are equivalent to small 
values of k for m = modes, the Doppler shift experienced 
by m = 2 modes means that rotational suppression instead 
occurs for k ~ 2f2*/fi or b, which is 84x (day/P*) for KOI-54's 
orbital period of 42 days. Fig. [7]b demonstrates this, where 
e.g. little power can be observed near k = 84 for P„ — 1 day. 

Furthermore, rotational suppression does not act on 

B It is sufficient to consider only nonnegative k, i.e. to use a uni- 
modal Fourier series, since the Fourier coefficient of orbital har- 
monic k must be the complex conjugate of that for — fc, since the 
lightcurve is real valued. 



low-fe harmonics in edge-on systems, as Fig.[7]b also shows. 
Indeed, since fast rotation Doppler shifts lower-order retro- 
grade modes — which radially couple more strongly to the 
tidal potential — toward values of k nearer to the Hansen 
peak of ~ 2f2 pcr i/f2 or b, the power in lower harmonics can 
even be enhanced by sufficiently fast rotation rates. 

The rightmost panel of Fig. [7]b shows the harmonic 
decomposition of our simple equilibrium tide model for an 
edge-on orientation, not including irradiation (§ 6.1 Ap- 



pendix B2 1 . Comparing this plot to the left four panels 



shows in particular that, in edge-on orbits, the dynamical 
tide is not rotationally suppressed for harmonics where the 
equilibrium tide has large amplitudes, unlike for face-on ori- 
entations. Thus the ellipsoidal variation of edge-on systems 
may be buried beneath the dynamical tidal response. This 
implies that full dynamical modeling may be necessary to 
constrain system parameters for edge-on binaries, and that 
care must be taken in searches for eccentric binaries, since it 
cannot be assumed that their lightcurves will be dominated 
by ellipsoidal modulations. 



6.4 Lightcurve power spectrum modeling 

We performed preliminary quantitative modeling of the pul- 
sation data in Table [2] As noted before, tidally driven pul- 
sations should have frequencies which are pure harmonics of 
the orbital frequency, ui — fcf2 or b for k £ Z. Although most 



© 0000 RAS, MNRAS 000, 000-000 



14 Burkart, Quataert, Arras, and Weinberg 



of the pulsations Wll report are of this form, some clearly 
are not, and are as such unaccounted for in linear perturba- 
tion theory. Hence we only attempted to model pulsations 
within 0.03 in k of a harmonic (set arbitrarily); this lim- 
ited our sample to 21 harmonics, as shown in Table [2] We 
provide an explanation for the nonharmonic pulsations in 

There are eight primary parameters entering into our 
modeling of the remaining observed harmonics: stellar 
masses Mi,2, radii Pi, 2, ZAMS metallicities Zi,2, and rigid- 
body rotation periods Pi, 2- To explore a range of stellar 



parameters, we used the stellar evolution code MESA (Pax- 



ton et al. 20111 to create two large sets of stellar models, 
one for each star, with ranges in M and R determined by 
Wll's constraints (Table [TJ. We set both stars' metallic- 
ities to 0.04. The other two parameters, Pi and P2, were 
treated within our nonadiabatic code using the traditional 
approximation. We set Pi = P2 = 1.5 days, comparable to 



the expected pseudosynchronous rotation period (§ 5.1 1 and 
qualitatively consistent with the small-amplitude flux per- 
turbations of lower harmonics seen in KOI-54 (§ 6.3 1. We 



fixed all of the orbital parameters to the mean values given 
in Wll. 

As discussed in § |6.6| it is possible that the 90th and 
91st harmonics observed in KOI-54 are m = ±2 modes re- 
sponsible for resonance locks, and are thus in states of nearly 
perfect resonance. Indeed, even if they are m = chance 
resonances, which are ~ 200 times easier to observe with 
K OI-5 4's face-on orbital inclination than m = ±2 modes 
(§ 3.3 1, we find that a detuning of |<5w/f2 rb| ~ 10 -2 is re- 
quired to reproduce the amplitude of either harmonic, where 
Suj is the difference between the eigenmode and driving fre- 
quencies (with 8uj — representing a perfect resonance). 

Such a close resonance represents a precise eigenfre- 
quency measurement, and should place stringent constraints 
on stellar parameters. However, this degree of resonance is 
also very difficult to capture in a grid of stellar models be- 
cause even changes in (say) mass of AM/M ~ 10~ 4 can alter 
the mode frequencies enough to significantly change the de- 
gree of resonance; future alternative modeling approaches 
may obviate this difficulty (§ [71. A second problem with 
trying to directly model the 90th and 91st harmonics is that 
the amplitudes of both of these harmonics may be set by 
nonlinear processes, as addressed in § |6.5| If correct, this 
implies that these particular modes strictly cannot be mod- 
eled using the linear methods we focus on in this paper. 

We are thus justified in restricting our analysis to only 
those integral harmonics in the range 35 ^ k ^ 89. We chose 
k = 35 as our lower bound to avoid modeling harmonics that 
contribute to ellipsoidal variation. We set m — for all of 
our analysis for the reason stated above. We also only used 
I — 2 for the tidal potential, since additional I terms are 
suppressed by a) further powers of P/P pGr i ~ 0.16, and b) 
smaller disk-integral factors from § 3.3 (e.g., 63/62 = 0.2). 

To find a reasonable fit to the harmonic power observed 
in KOI-54, we attempted a simplistic, brute-force optimiza- 
tion of our model against the data: we first modeled the lin- 
ear response of each stellar model in our grid separately, ig- 
noring its companion, and calculated the resulting observed 
flux perturbations as a function of k. We then compared 
the absolute values of these flux perturbations to the ob- 
servations of KOI-54 and selected the best 10 3 parameter 
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Figure 8. Exploratory modeling of the Fourier power spectrum 
of KOI-54's lightcurve. Our inhomogeneous theoretical model, 
including nonadiabaticity and rotation in the traditional approx- 
imation, is plotted above the graph's horizontal axis, while the 
data for KOI-54 (Table ^ is plotted below. Parameters cor- 
responding to this plot are in Table [4] Since harmonics 90 & 
91 represent extreme resonances, they are difficult to resolve in 
a given stellar model grid, and fits which reproduce their am- 
plitudes cannot reproduce other parts of the Fourier spectrum. 
Thus we attempted to fit only 35 ^ k ^ 89 and not the shaded 
region. Our fitting process was simplistic (see text), and we did 
not approach a full optimization, although our best fit does agree 
reasonably well with the data. Fig. [9] shows the observed flux 
perturbation from this plot for both stars separately on a fine 
frequency grid. 



Table 4. Stellar parameters used in Figures [8] and [9] 



star 


M/Mq 


R/Rq 


Z 


P»/day 


1 


2.278 


2.204 


0.04 


1.5 


2 


2.329 


2.395 


0.04 


1.5 



choices (Af , R) for each star. (In future work, pulsation 
phases should be modeled in addition to the amplitudes re- 
ported by Wll, since this doubles the information content 
of the data; see also further discussion of phases in 



4.2 



) 

Given this restricted set of stellar models, we computed the 
theoretical Fourier spectra for all 10 6 possible pairings of 
models. 

Fig. [8] shows one of our best fits to the observations of 
KOI-54; Table [1] gives the associated stellar parameters. We 
obtained many reasonable fits similar to Fig. [8] with dissim- 
ilar stellar parameters, demonstrating that many local min- 
ima exist in this optimization problem. As a result, Fig. [8] 
and Table [4] should not be interpreted as true best fits but 
rather as an example of a model that can semi-quantitatively 
explain the observed harmonic power in KOI-54. We leave 
the task of using the observed pulsation data to quantita- 
tively constrain the structure of the stars in KOI-54 to future 
work, as we discuss in § [7] 

Responses from both stars were used to create the plot 
in Fig. [8] Fig. [9] on the other hand, uses the same param- 
eters, but instead shows each star's observed flux pertur- 
bation separately and evaluated on a fine grid in frequency 
rather than only at integral orbital harmonics. As a result, 
Fig. [9] exposes the position of normal modes (which corre- 
spond to peaks in the black curves) in relation to observed 
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Figure 9. Individual stars' contributions to the Fourier spectrum from Fig. [8] (black curves), evaluated on a fine grid in k = u/Q OI yj. 
The actual Fourier series decomposition (Fig.|8j is obtained by adding both stellar responses (together with phases) at each integer value 
of k. Observed harmonics are displayed as red vertical lines, with short red horizontal lines indicating their amplitudes. We used a mesh 
of 50 points per unity increment in k to produce these plots, and find that the highest peaks of the eigenmode resonances nearest k ~ 90 
are then at the same amplitude level as the observed 90th and 91st harmonics; this means that a detuning of Ak = Slu/Q^^, ~ 0.02 is 
required to explain these pulsations if they have m = (see text). The regions below the minimum amplitude reported by Wll (Table[2]l 
are shaded. 



harmonics (shown as red vertical lines), as well as other fea- 
tures not captured in Fig. [8]s raw spectrum. Fig. [9] also 
shows that harmonics 90 and 91 must come from different 
stars if they are indeed m = g-modes (although this may 
not be the case; see § |6.6[ l , since the g-mode spacing near 
k ~ 90 is much larger than the orbital period (with the same 
logic applying for harmonics 71 and 72). 



6.5 Nonharmonic pulsations: three-mode coupling 

Wll report nine pulsations which are not obvious harmonics 
of the orbital frequency; these have asterisks next to them in 
Table[2] As we showed previously (§ 4.2 and Appendix Al I, 
these cannot be linearly driven modes. Here we present one 
possible explanation for the excitation of these pulsations. 

To begin, we point out the following curious fact: the 
two highest-amplitude nonharmonic pulsations in Table [2] 
(F5 and F6) have frequencies which sum to 91.00 in units 
of the orbital frequency — precisely the harmonic with the 
second-largest amplitude (F2). (This is the only such in- 
stance, as we discuss below.) 

Although this occurrence could be a numerical coinci- 
dence, it is strongly suggestive of parametric decay by non- 
linear three-mode coupling, the essential features of which 
we now describe. First, however, we emphasize that the 
treatment we present here is only approximate. In real- 
ity, the process of nonlinear saturation is much more com- 
plicated, and a more complete calculation would involve 
fully coupling a large number of eigenmodes simultaneously 
( Weinberg fc Quataert|2008 1. 

If a parent eigenmode is linearly excited by the tidal 
potential to an amplitude that surpasses its three-mode- 
coupling threshold amplitude S a , any energy fed into it 
above that value will be bled away into daughter mode pairs, 
each with frequencies that sum to the parent's oscillation fre- 
quency ( Weinberg et al.||2011 l. In a tidally driven system, 
the sum of the daughter modes' frequencies must thus be a 
harmonic of the orbital frequency. 



For a parent with indices a = (n,l,m) linearly driven 
at a frequency er, the threshold is given by 



where 



I Sal 



mm 

be 



{Tabcj , 



ALUbOJc\Kabc 



(37) 



(38) 



uji is a mode frequency, 7^ is a mode damping rate, n a bc is 
the normalization-dependent nonlinear coupling coefficient 

, is the detuning 



(Schenk et al. 20021, 5uj b 



U>b 



frequency, and the minimization is over all possible daughter 
eigenmodes b and c (each short for an (n, I, m) triplet) (j The 
nonlinear coupling coefficient K a bc is nonzero only when the 
selection rules 



= mod(/ a + l b + l c , 2), 
= m a + m b + m c , 

lib — lc < la < lb + lc 



(39) 
(40) 
(41) 



are satisfied. Due to the second of these rules, any Doppler 
shifts due to rotation do not affect the detuning since they 
must cancel. 

For a simple system of three modes, the nonlinear cou- 
pling's saturation can be determined analytically. The par- 
ent saturates at the threshold amplitude S a , and the ratio 
of daughter energies within each pair is given by the ratio 
of the daughters' quality factors: 

Eb _ qb _ Ub/jb , 42 s 
E c q c u c /j c 

Equations ( |37[ ) and ( |38[ ) exhibit a competition that 
determines which daughter pair will allow for the lowest 
threshold. At larger daughter /, modes are more finely 



7 This section uses the normalization of |Weinberg et al.| |201l[|, 
whereas the rest of the paper uses the normalization given in § |3.2| 
We of course account for this when giving observable quantities. 
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spaced in frequency, since g-mode frequencies roughly sat- 
isfy the asymptotic scaling from eq. hence, the detuning 
Subc becomes smaller (statistically) with increasing I. How- 
ever, higher daughter I also leads to increased damping rates 
at fixed frequency 



17) 



As such, the minimum threshold 
will occur at a balance between these two effects. 

In order to semi-quantitatively address the phenomenon 
of three-mode coupling in KOI-54, we produced an example 
calculation of S a together with a list of best-coupled daugh- 
ter pairs. To this end, we used a MESA stellar model ( |Pax-| 
ton et aL||2011 \ consistent with the mean values of star l's 
properties reported in Wll (Table [T]). We computed this 
model's adiabatic normal modes using the ADIPLS code 



( Christensen-Dalsgaard 20081, and calculated each mode's 



global quasiadiabatic damping rate 7„; due to radiative dif- 
fusion (§[3^2|. 

We focus on the second-highest-amplitude k — 91 
harmonic present in the data (Fl from Table [2j and set 
a = 91 x il or b; as pointed out in Wll, for m a = the 
quadrupolar eigenmode with natural frequency closest to 
the 91st orbital harmonic is the gi4 mode, i.e., the g-mode 
with 14 radial nodes. We thus take this as our parent mode. 



The minimization in eq. ( 37 1 is over all normal modes, 



of which there is an infinite number. To make this problem 
tractable numerically, we essentially followed the procedure 



described in |Weinberg fc Quataert ( 2008 1 



(i) We restricted daughter modes to 1 ^ I ^ 6. There 
is no reason a priori to suggest I should be in this range, 
but, as shown in Table [5] 1 ^ I ^ 3 turns out to be the 
optimum range for minimization in this particular situation, 
and modes with I > 6 are irrelevant. 

(ii) The quantity to be minimized in eq. (37 1, T a bc, 
achieves its minimum at fixed Sujb c and n a bc for u>b ~ u c , 



given the scaling from eq. (17 1. As such, we computed all 



normal modes b with frequencies in the range / < uJb/ijJ a < 
1 - /; we took / = 1/10, which yielded 344,479 potential 
pairs, but trying / = 1/5, which yielded 61,623, did not 
change the result. 

(iii) We computed T a bc, not including the three-mode- 
coupling coefficient K a b c (since it is computationally expen- 
sive to evaluate), for all possible pairs of modes satisfying 



(i) and (ii) as well as the selection rules in equations ( 39 1 — 
(4T|. 

(iv) From the results of (iii), we selected the N = 5000 
smallest threshold energies, and then recomputed T a bc for 
these pairs this time including K a bc ( Weinberg et al.||2011 |. 
(Trying N = 1000 did not change the results.) We set mb = 
m c = = m a for simplicity, since n a bc depends only weakly 
on the values of m so long as eq. ( 40 1 is satisfied. Sorting 



again then yielded the best-coupled daughter pairs and an 
approximation for the saturation amplitude S a - 

Table [5] shows the best-coupled daughter mode pairs 
resulting from this procedure. It is interesting to note that 
most daughter pairs a) involve an I = 1 mode coupled to an 
I = 3 mode (P2, P4, P5, P6), and/or b) have a large quality- 
factor ratio (all except P2, P5, & P9 have || log 10 (qb/ q c )\ > 
0.5). 

For daughter pairs satisfying (a), the I = 3 mode would 
be much harder to observe in a lightcurve since disk averag- 
ing involves strong cancellation for larger-Z modes — indeed, 
Table [3] shows 63/61 ~ 0.1 for Eddington limb darkening, 



where bi is a disk-integral factor defined in eq. (231. (The 
other disk-integral factor, ci, does not decline as sharply 
with increasing I, but corresponds to cross-section perturba- 
tions, which are small relative to emitted flux perturbations 
as discussed below.) For daughter pairs satisfying (b), since 
the ratio of daughter amplitudes scales as the square root 
of the ratio of their quality factors, one of the modes would 
again be difficult to observe. 

Furthermore, if the parent had m a = ±2 instead of 
m a = (see § 6.6 1, each daughter pair would have sev- 
eral options for rrib and m c , introducing the possibility 
of \rrib\ 7^ \m c \. This would mean daughters would ex- 
perience even greater disparity in disk-integral cancella- 
tion due to the presence of Yi r 

\Y 10 (e o ,4> o )/Y 32 (e a ,<i> o )\ 



1221; e.g., 



<f>o) m eq. 
0.02 for KOI-54. 
The above results provide a reasonable explanation for 
why there is only one instance of two nonharmonic pulsa- 
tions adding up to an observed harmonic in the data for 
KOI-54 — only P9 from Table [5] has the potential to mimic 
pulsations F5 and F6 from Wll. Nonetheless, the nonlin- 
ear interpretation of the nonharmonic pulsations in KOI-54 
predicts that every nonharmonic pulsation should be paired 
with a lower-amplitude sister such that their two frequen- 
cies sum to an exact harmonic of the orbital frequency. This 
prediction may be testable given a sufficient signal-to-noise 
ratio, which may be possible with further observations of 
KOI-54. 

Lastly, we can attempt to translate our estimate of the 
parent threshold amplitude S a into an observed flux per- 
turbation, <5J sa t/J, using the techniques of § 3.3 Since our 
nonlinear saturation calculation was performed with adia- 
batic normal modes, we strictly can only calculate the ob- 
served flux variation due to cross-section perturbations, 8 J cs 



(the £ r component of eq. 22 1, and not that due to emitted 
flux perturbations, 8J a { (the AF component of eq. 221. It 
evaluates to 

5 J cs 



\S a x (26/ -c,) x Zr,a(R) x Y 20 (e o ,(j> o )\ 
~ 1.7 mmag, 

However, we can employ our nonadiabatic code to calibrate 
the ratio of 5J CS to 8J e i , which we find to be SJ e f/SJ cs — 9 for 
the 91st harmonic. We can then estimate the total saturated 
flux perturbation: 

/if. \ XT 

17 mmag. 



Jsat 


- ( SJ - ! 


►0 


5 J cs 


J 


\3Jcb 


J 



This result is a factor of ~ 100 too large relative to 
the observed amplitude of 229 ^imag for the 91st harmonic 
(Table [2]|. Taken at face value, this would mean that the 
inferred mode amplitude is below threshold, and should not 
be subject to nonlinear processes, despite evidence to the 
contrary. There are several possible explanations for this 
discrepancy. If the 91st harmonic is actually an m — ±2 
mode, which we proposed in § |5.1[ then the intrinsic ampli- 
tude required to produce a given observed flux perturbation 
is a factor of ~ 200 times larger than for m = modes given 
KOI-54's face-on inclination (§ |3.3[ ). This would make the 
observed flux perturbation of the 91st harmonic compara- 
ble to that corresponding to the threshold for three-mode 
coupling, consistent with the existence of nonharmonic pul- 
sations in the lightcurve. We discuss this further in § |6.6| 
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Table 5. Ten best-coupled daughter mode pairs resulting from the procedure outlined in steps (i - iv) of § |6,5| This is an example 
calculation and is not meant to quantitatively predict the nonharmonic components of KOI-54's lightcurve. All frequencies and damping 
rates are in units of £1 or t,. The square root of the daughter quality factor ratio, \J Qb/qc, gives an estimate of the ratio of daughter mode 
amplitudes, and hence of their potential relative lightcurve contributions. 



ID 


(h 


,n b ) 


(lc 


n c ) 


^6 




§ l°glo(<?b/<7c) 


log 10 \Suj bc \ 


l°gio( 2 v 


PI 


(2, 


-37) 


(2, 


-23) 


35.3 


55.8 


0.66 


-1.4 


-2.1 


P2 


(1. 


-25) 


(3, 


-30) 


29.9 


61.0 


0.057 


-1.5 


-2.5 


P3 


(1. 


-28) 


(1. 


-11) 


26.8 


64.1 


0.69 


-1.3 


-3.1 


P4 


(1. 


-50) 


(3, 


-24) 


15.3 


75.7 


1.2 


-1.4 


-1.7 


P5 


(1. 


-27) 


(3, 


-29) 


27.8 


63.2 


0.031 


-1.4 


-2.5 


P6 


(1. 


-42) 


(3, 


-25) 


18.0 


73.1 


1.2 


-1.6 


-1.7 


P7 


(1, 


-36) 


(1, 


-10) 


20.9 


69.8 


1.3 


-0.61 


-2.7 


P8 


(2, 


-35) 


(2, 


-24) 


37.2 


53.7 


0.51 


-0.87 


-2.2 


P9 


(2, 


-29) 


(2, 


-28) 


44.7 


46.4 


0.038 


-0.99 


-2.5 


P10 


(1. 


-35) 


(1. 


-10) 


21.5 


69.8 


1.2 


-0.51 


-2.8 



Alternatively, if the 91st harmonic is in fact an m = 
mode, many daughter modes may coherently contribute to 
the parametric resonance, reducing the threshold consider- 
ably, as in Weinberg et al. (20111. A more detailed cal- 



culation, coupling many relevant daughter and potentially 
granddaughter pairs simultaneously, should be able to ad- 
dress this more quantitatively. 



6.6 Are harmonics 90 and 91 caused by prograde, 
resonance- locking, \m\ — 2 g- modes? 

As introduced in § |5.1| having two pseudosynchronized stars 
presents an ostensibly appealing explanation for the large- 
amplitude 90th and 91st harmonics observed in KOI-54 
(henceforth Fl and F2; Table|2|: each is the manifestation of 
a different highly resonant eigenmode effecting a resonance 
lock for its respective star by opposing the equilibrium tide's 
torque. 

We discuss the viability of this interpretation below. 
First, however, what alternate explanation is available? The 
most plausible would be that Fl and F2 are completely inde- 
pendent, resonantly excited m — modes. Each coincidence 
would require a detuning of \uj n i — o"fc m |/f2orb ~ 2 x 10~ 2 
(Fig. [jjj, which is equivalent to \oj n i — <7fc TO |/w n i ~ 10 -4 , 
where u n i is the nearest eigenfrequency and G^ m = fcf2 or b — 
mfi* is the driving frequency. The probability of having a 
detuning equal to or smaller than this value, given ~ 10 
available modes (Fig. |9j, is ~ 10%, so the combined proba- 
bility if the resonances are independent is ~ 1%. Moreover, 
in § |6.4| we show that in this m — interpretation, Fl and F2 
must come from different stars, yet there is no explanation 
for why the two excited modes are so similar. 

If instead Fl and F2 are due to highly resonant m = ±2 
resonance locking modes, several observations are naturally 
explained. The high degree of resonance is an essential fea- 
ture of the inevitable pseudosynchronous state reached when 
the torque due to the dynamical tide cancels that due to the 
equilibrium tide (§ |5.1[ ). The fact that the resonant modes 
correspond to similar k would be largely a consequence of 
the fact that the two stars in the KOF54 system are sim- 
ilar in mass and radius to ~ 10%, so that a similar mode 
produces the dynamical tide torque in each star (although a 



corresponding ~ 10% difference in k would be equally pos- 
sible in this interpretation). 

In addition, we showed in § |6.5| that the observed ampli- 
tudes of Fl and F2 are a factor of up to ~ 100 smaller than 
their nonlinear threshold values assuming m = 0. There is 
also strong evidence that at least F2 has its amplitude set by 
nonlinear saturation. Having m 7^ would help to resolve 
this discrepancy because the intrinsic amplitude of m = ±2 
modes would need to be ~ 200 times larger to produce the 
observed flux perturbation. This would then imply that the 
amplitudes of Fl and F2 are indeed above the threshold for 
three-mode coupling, naturally explaining the presence of 
the nonharmonic pulsations in the KOI-54 lightcurve. 

However, several significant problems with the 
resonance-locking interpretation arise upon closer exami- 
nation. Assume that Fl and F2 indeed correspond to 
m = ±2 g-modes that generate large torques effecting 
P ps ~ 1.8 days pseudosynchronization locks. In order to 



create positive torques, eq. |C4| shows that we must have 
m(kQ or b — mS7«) > 0, which reduces to 



(fc/m)n or b > fi*. 



(43) 



In order to determine which modes correspond to Fl and 
F2, we can enforce a close resonance by setting 



: 90 n„ r 



2Q, 



(44) 



where we have used the fact that eq. ( |43[ ) requires k and m to 
have the same sign for a positive torque. For / = 2 and using 
a MESA stellar model consistent with Wll's mean modeled 
parameters for star 1 (Table [l]), eq. (44 1 yields n ~ 30, 
neglecting rotational modification of the modes (i.e., not 
employing the traditional approximation). 

However, in our calculations in § |5.1| we find that the 
resonant torque due to the dynamical tide is instead typi- 
cally caused by g-modes with n of 8 - 15 (basically set by 
the intersection of the Hansen coefficient and linear overlap 
curves, as discussed in § |4.2| in the context of flux pertur- 
bations). Using eq. (44 1 again, this would mean we would 
expect k of 140 - 200. Furthermore, we find that even a 
perfectly resonant n = 30 g-mode makes a negligible contri- 
bution to the torque. This is true both for ZAMS models 
and for evolved models consistent with the observed radii in 
KOI-54, indicating that there is little uncertainty introduced 
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by the details of the stellar model. This result suggests that 
the g-modes inferred to correspond to Fl and F2 are in- 
consistent with what would be expected from our torque 
calculation if the rotation rate is indeed ~ 1.8 days. 



If we account for rotation in the traditional approxima- 



tion 



3.4 1, the n of a prograde mode of a given frequency 
4/6 times smaller than its cor- 



can be at most a factor of 
responding nonrotating value; this follows from the fact that 
the angular eigenvalue A asymptotes to m 2 = 4 in the limit 
uj n i <C fi* for prograde modes, instead of A = 1(1 + 1) = 6 
in the nonrotating limit. This reduces the n of the 90th 
harmonic from n ~ 30 to n ~ 24, still insufficient to yield a 
significant torque. 

Another major problem with the resonance lock inter- 
pretation is that although our orbital evolution simulations 
described in § |5 . l| ubiquitously produce resonance locks, they 
always occur in only one star at a time. This is because if 
a mode is in a resonance lock in one star and a mode in the 
other star tries to simultaneously resonance lock, the first 
lock typically breaks since the orbital frequency begins to 
evolve too quickly for the lock to persist. Although it is 
possible for simultaneous resonance locks in both stars to 
occur, such a state is very improbable. Similar orbital evo- 



lution simulations presented in Fuller & Lai (2011 1 did pro- 



duce simultaneous resonance locks, but only because they 
simulated only one star and simply doubled the energy de- 
position rate and torque, thus not allowing for the effect just 
described. 



Finally, we point out one last inconsistency in the res- 
onance lock interpretation of Fl and F2. It is straightfor- 
ward to calculate the predicted flux perturbation associated 
with perfectly resonant \m\ = 2 g-modes in resonance locks 
(using, e.g., the calibration discussed at the end of § 6.51: 
for modes ranging from n ~ 8 - 15, we find that the pre- 
dicted flux perturbation for KOI-54's parameters is ~ 10 - 
30 ^imag. This is a factor of ~ 10 smaller than the observed 



flux perturbations, yet somewhat larger than the smallest- 
amplitude pulsation reported by Wll. It is also a factor 
of ~ 2 smaller than the nonlinear coupling threshold for an 
m — 2 mode (which we determined using the same proce- 
dure as in § 6.5 extended to allow for an m 7^ parent), 



although the uncertainties involved in our nonlinear esti- 
mates are significant enough that we do not consider this to 
be a substantial problem. 

Thus even if Fl and F2 can be attributed to modes un- 
dergoing resonance locks (which is highly nontrivial, as we 
have seen), the observed amplitudes are larger than those we 
predict. Conversely, if Fl and F2 are simply chance m — 
resonances, it appears that if a resonance lock existed, it 
would have been detected, although the possibility exists 
that the resonant mode's flux perturbation was marginally 
smaller than those of the 30 reported pulsations due to un- 
certainties in our calculations. Firmer constraints on the 
flux perturbations in KOI-54 at k ~ 140 - 200 would be 
very valuable in constraining the existence of such m = ±2 
modes, as would information about the phases of the 90th 



7 DISCUSSION 

We have developed a set of theoretical tools for understand- 
ing and modeling photometric observations of eccentric stel- 
lar binaries. This work is motivated by the phenomenal 
photometry of the Kepler satellite and, in particular, by the 
discovery of the remarkable eccentric binary system KOI-54 
I Welsh et al.||2011 henceforth Wll). This system consists 
of two similar A stars exhibiting strong ellipsoidal lightcurve 
variation near periastron passage due to the system's large 
(e = 0.83) eccentricity. Wll successfully modeled this phe- 
nomenon, and also reported the detection of at least 30 dis- 
tinct sinusoidal pulsations in KOI-54's lightcurve (§[5|, ~ 20 
at exact harmonics of the orbital frequency and another ~ 
10 nonharmonic pulsations. Although our work has focused 
on modeling KOI-54, our methods and techniques are more 
general, and are applicable to other similar systems. 

We developed a simple model of KOI-54's periastron 
brightening, including both the irradiation and equilibrium 
tide components of this effect, which agrees at the ~ 20% 
level with the results Wll obtained using a much more de- 



tailed simulation (§ 6.1 1. Our model may be useful for analy- 



sis of other eccentric stellar binaries, allowing determination 
of orbital and stellar parameters; its simplicity should enable 
it to be implemented in an automated search of Kepler data. 

In § |4| we used the adiabatic normal mode formalism 
(see § |3.2| and, e.g., Christensen-Dalsgaard 2003 Kumar 



et al. 19951, to establish a qualitative connection between 



the range of stellar modes excited in a given binary system 
and the system's orbital properties. For more detailed quan- 
titative modeling of the harmonic pulsation spectrum of a 
given binary system, we further developed the nonadiabatic, 



inhomogeneous tidal method from Pfahl et al. ( 2008 1 by in- 



cluding the Coriolis force in the traditional approximation 
(§|6.2| Appendix |A). 

In § |6.3| we used this method to show that fast rota- 
tion tends to suppress power in the lower harmonics of a 
face-on binary system's lightcurve (Fig. FFJ. This can qual- 
itatively explain why there is a scarcity of large-amplitude, 
lower-harmonic pulsations in KOI-54's lightcurve, relative to 
predictions for nonrotating stars (Fig. [3J . We also showed 
in § |6.3[ however, that the dynamical tidal response may be 
much larger than ellipsoidal variation in edge-on binaries, 
unlike in KOI-54 (which has an inclination of i — 5.5°; see 
Table [I]). For such systems, simultaneous modeling of the 
dynamical and equilibrium tides may be required in order 
to constrain system properties. 

Moreover, in §[5] we showed that rapid rotation periods 
of ~ 1.8 days are expected for the A stars in KOI-54, due 
to pseudosynchronization with the orbital motion near pe- 
riastron. This pseudosynchronous rotation period is shorter 
than the value of 2.53 days assumed by Wll. The latter 
value is appropriate if the only appreciable torque is that 



and 91st harmonics (see § 4.2 1 



produced by the equilibrium tide (Appendix C2 1 . Since res- 
onantly excited stellar g-modes can produce a torque com- 
parable to that of the equilibrium tide, pseudosynchronous 
rotation can occur at even shorter rotation periods (Fig. B. 
This involves a stochastic equilibrium between prograde res- 
onance locks and the equilibrium tide. These same rapid ro- 
tation periods (~ 1.8 days) yield predicted lightcurve power 
spectra that are the most qualitatively consistent with the 
pulsation data for KOI-54 (Fig. [7]). 
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In § |6.4| we performed a preliminary optimization of our 
nonadiabatic model by comparing its results in detail to the 
Fourier decomposition of KOI-54's lightcurve (Table [2|. We 
searched over an extensive grid of stellar masses and radii, 
assuming a metallicity of twice solar and a rotation period 
of 1.5 days. We also set m = 0, since KOI-54's nearly face- 
on orientation implies that this is the case for almost all of 



the pulsations we modeled (§ 3.3 1. The modeling challenge 



in tidal asteroseismology contrasts with that of standard 
asteroseismology in that a) we must simultaneously model 
both stars, and b) pulsation amplitudes and phases contain 
the key information in our case, since we are considering 
a forced system, whereas pulsation frequencies constitute 
the data in traditional asteroseismology. Moreover, stellar 
rotation is sufficiently rapid in eccentric binaries that its 
effect on stellar g-modes cannot be treated perturbatively. 

Although our minimization procedure was quite sim- 
ple, we were able to obtain stellar models with power spec- 
tra semi-quantitatively consistent with the observations of 
KOI-54 (Fig. [8]& Fig. [9). The resulting model in Fig. [|] is 
not formally a good fit, but this is not surprising given that 
two of the key parameters (metallicity and rotation period) 
were not varied in our analysis. Moreover, in our preliminary 
optimization we found that there were many local minima 
that produced comparably good lightcurves. 

As noted above, a priori calculations suggest that both 
stars in the KOI-54 system should have achieved a pseu- 
dosynchronous state at rotation periods of ~ 1.8 days. This 
requires frequent resonance locks to occur, when a single 
\m\ = 2 eigenmode comes into a near-perfect prograde reso- 
nance. A natural question is whether such a highly resonant 
mode could contribute to the KOI-54 lightcurve; this possi- 
bility is particularly attractive for the two largest-amplitude 
harmonics observed, the 90th and 91st. (See also our calcu- 
lation of nonlinear saturation from 



6.5 discussed below. 



However, we find quantitative problems with this inter- 
pretation (§ 6.6 1. First, our orbital evolution simulations 



(§ 5.1 1 indicate that only one resonance lock should exist at 



a time, meaning that only one of the two large-amplitude 
harmonics could be explained in this way. This result is in 
disagreement with the simulations performed by |Fuller fe| 
|Lai| p0lT| >, since they did not simultaneously model both 
stars. 

Further, in our calculations, the g-modes capable of pro- 
ducing torques large enough to effect resonance locks have 
n typically in the range of 8 - 15 (where n is the number 
of radial nodes), while the 90th harmonic corresponds to n 
of 25 - 40 for m = ±2 and rotation periods of 2.0 - 1.5 
days. Also, we predict that g-modes producing resonance 
locks should have k of 140 - 200, much larger than ~ 90, 
and flux perturbations of 10 - 30 ^tmag. The latter values 
are a factor of ~ 10 less than that observed for the 90th 
and 91st harmonics, but slightly larger than the smallest 
observed pulsations. 

It thus seems quantitatively difficult to interpret har- 
monics 90 and 91 in KOI-54 as manifestations of m = ±2 
modes in resonance locks, although we cannot conclusively 
rule out this possibility. Instead, it seems likely that they 
are simply chance m = resonances (as is almost certainly 
the case for the overwhelming majority of the other observed 
pulsations in KOI-54). One theoretical uncertainty resides 
in our omission of rotational modification of the stellar eigen- 



modes when computing tidal torques. Our estimates suggest 
that this is a modest effect and is unlikely to qualitatively 
change our conclusions, but more detailed calculations are 
clearly warranted. 

We note that in future work, pulsation phases should 
be modeled in addition to the amplitudes reported by Wll, 
since this effectively doubles the information content of the 
data. Indeed, we showed in § |4.2| that a resonant pulsation's 
phase is strongly influenced by the mode's value of m. In 
particular, since harmonics 90 and 91 are likely standing 
waves, as can be seen in the propagation diagram in Fig. [I] 
measurement of their phases could help to resolve the uncer- 
tainties pointed out above by supplying direct information 
about their degrees of resonance, thus potentially confirming 
or disproving the m = ±2 resonance lock interpretation. 

In § |6.5| we pointed out evidence for nonlinear mode 
coupling in KOI-54's observed pulsations: the existence of 
nonharmonic pulsations (which does not accord with lin- 
ear theory; § 4.2 1 and the fact that two of them have fre- 



quencies that sum to exactly the frequency of the 91st har- 
monic, the second-largest-amplitude harmonic pulsation in 
KOI-54's lightcurve. This is consistent with parametric reso- 
nance, the leading-order nonlinear correction to linear stellar 
oscillation theory (Weinberg et al.|2011 |. 

Motivated by this observation, we performed a nonlin- 
ear stability calculation that qualitatively explains why no 
other similar instance of a nonharmonic pair summing to an 
observed harmonic is present in the data: for the majority 
of daughter pairs likely to be nonlinearly excited, there are 
sufficient differences in the I and m values of the daughter 
pair members, or sufficient differences in their predicted sat- 
urated energies, that only one member of the pair would be 
observable given current sensitivity. Nonetheless, the non- 
linear interpretation makes the strong prediction that ev- 
ery nonharmonic pulsation should be paired with a lower- 
amplitude sister such that their two frequencies sum to an 
exact harmonic. This prediction may well be testable given 
a better signal-to-noise ratio. 

One additional feature of the nonlinear interpretation 
is that if the nonlinearly unstable parent is an m = mode, 
then the threshold amplitude for a linearly excited mode to 
be unstable to parametric resonance, which we have just ar- 
gued exists in KOI-54, implies flux perturbations that are 
a factor of ~ 100 larger than those observed. In contrast, 
the parent being an m = ±2 mode ameliorates this discrep- 
ancy because the parent's intrinsic amplitude must be ~ 200 
times larger for a given flux perturbation due to KOI-54's 
face-on orientation. This result thus argues in favor of the 
91st harmonic in KOI-54 being an m — ±2 mode caught in 
a resonance lock, as discussed above. 

There are many prospects for further development of 
the analysis begun in this paper. For example, in traditional 
asteroseismology, standard methods have been developed al- 
lowing a set of observed frequencies to be inverted uniquely, 
yielding direct constraints on stellar parameters, including 
the internal sound speed profile (Unno et al. 1989). The 



essential modeling difficulty in tidal asteroseismology is our 
inability to assign each observed pulsation amplitude to ei- 
ther star of a given binary a prion, hindering our attempts 
to develop a direct inversion technique. We leave the exis- 
tence of such a technique as an open question. 

Future observations of eccentric binaries may avoid this 
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difficulty if one star is substantially more luminous than the 
other. However, for eccentric binaries with similar stars, 
in the absence of a means of direct inversion, we are left 
with a large parameter space over which to optimize, con- 
sisting at minimum of eight quantities: both stars' masses, 
radii/ages, metallicities, and rotation periods. Even this pa- 
rameter set may ultimately prove insufficient, if modeling of 
tidally forced pulsations is found to be sensitive to the de- 
tails of e.g. chemical mixing or convective overshoot, which 
can modify the Brunt- Vaisala frequency and thus g-mode 
frequencies. 

One possible approach that should be explored in fu- 
ture work is to apply standard numerical optimization algo- 
rithms such as simulated annealing to this parameter space, 
attempting to minimize the \ 2 °f our nonadiabatic code's 
theoretical Fourier spectrum against the observed harmonic 
pulsation data. In practice, it may be preferable to develop 
interpolation techniques over a grid of models given the high 
resolution in stellar parameters needed to resolve the close 
resonances responsible for large-amplitude pulsations. 

Although KOI-54's stars lie near the instability strip, 
this fact is unimportant for the tidal asteroseismology theory 
presented in this work. Consequently, future high-precision 
photometric observations of other eccentric binaries may 
supply a window into the structure of stars previously inac- 
cessible by the techniques of asteroseismology. Constructing 
a data pipeline capable of reliably flagging eccentric binary 
candidates — e.g., finding efficient ways of searching for the 
equilibrium tide/irradiation lightcurve morphologies shown 
in Fig. |5| (Appendix [B| — is also an important, complemen- 
tary prospect for future work. 
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APPENDIX A: NONADIABATIC TIDALLY DRIVEN OSCILLATION EQUATIONS 



Here we will describe the computational procedure we employed to solve for tidally driven stellar responses, which we then 
used to model KOI-54's lightcurve. In Appendix |Al[ we account for rotation only by using Doppler-shifted driving frequencies 
fef2 or b — mQ* , and neglect any effects of the Coriolis force; in Appendix A2 we invoke the traditional approximation ( Bildsten 
et al.|1996) to account for the Coriolis force 



3.41. 



Al Formalism without the Coriolis force 

The gravitational potential due to the secondary, experienced by the primary, is given by 

GM 2 



U 2 ^i = - 



D — r\ 



(Al) 



Performing a multipole expansion ( Jackson|1999 \ and excising the I = (since it is constant) and I = 1 (since it is responsible 
for the Keplerian center-of-mass motion) terms, we are left with the tidal potential: 



v ' 1=2 m=-l 



\D(t) 



(A2) 



where 



4ir 
21 + 1 



yr™(7r/2,o) 



mod(Z +m + 1,2) 



/ 4tt 


(l + m- 


1)!! 


(l-m- 1)!! 


21 + 1 


(l + m 


!! 


{l-m)\\ 



(A3) 



Next, we shift to the primary's corotating frame (by sending 
terms of the Hansen coefficients: 



Mi 



i+i 



> <f> + Q t t) and expand the time dependence of the orbit in 
^2 WimY lm {e, 4>) eM-i<7kmt)Xf m (e)Ui(r), (A4) 



with Ohm = fcfiorb — mQ,* and 



Ut(r) 



( GMi 



(A5) 



The unit-normalized Hansen coefficients X lm were defined in eq. (Ill; here we are using the conventionally normalized Hansen 
coefficients Xf m = X^ m /(1 — e) l+1 , which are convenient to evaluate numerically as an integral over the eccentric anomaly: 

dE. 



Xi m = — i (1 — e cos E) cos 

7T ' 



k(E — e sin E) — 2m arctan 



1 + e 



1-e 



tan(£/2) 



(A6) 



If we represent the linear response of a star to the perturbing tidal potential by an abstract vector y(r,9,<f>,t) whose 
components are the various oscillation variables (e.g., (, r /r), then y can also be expanded, again in the primary's corotating 
frame, as in ( A4 1 : 



y = 



Ma 
M 1 



W l™ Y i™(0, <t>) Y ex P(- icr k m t)Xj c m (e)y^ n 



(A7) 



The equations necessary to determine y\ m (r ) are given in the appendix of Pfahl et al. ( 2008 \ , along with appropriate boundary 
conditions; note that their U is our Ui and their driving frequency to is our <Tfe m . 

After determining y\ m (r) in the corotating frame, we can switch to the inertial frame specified in 



3.1 



v ~ m7^ 



exp(-ikn olb t) W lm Yi m (6, 4>)Xf m (e)yf m (r 



(A8) 



As noted in § |4.2| we see in eq. \A&\ that the observed frequencies should be pure harmonics of the orbital frequency, even 
though the corresponding amplitudes of observed pulsations are influenced by the star's rotation rate (via the Doppler-shifted 
frequency cr fcm ). 



A2 Rotation in the traditional approximation 

We now invoke the traditional approximation (§ |3.4[ |; we must correspondingly adopt the Cowling approximation and employ 
the Hough functions (§ |5| as angular basis functions instead of spherical harmonics. 
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We expand the Hough functions as ( Longuet-Higgins 1968) 



= 2tt 



(A9) 



where P; m is a normalized associated Legendre function defined by 

Pim -■ 



2l + l(l-m) 



Y lm (O,0) = e^P lm (cos8) 



47r (l + m)\ 



(A10) 



We used the numerical method of calculating the expansion coefficients e\ lm detailed in Ogilvie & Lin (20041 § 5.4. The tidal 
potential in the corotating frame is then 



Mi 4^ 



(+1 



^ W; m e lm0 ^ exp(-i<7 fcm t)Xf m ^ H% m (n)Ux lm (r 



(All) 



where 



V fli / \Ri 



k 



<5km = fcf2 or b — mfi, , and the Coriolis parameter q on which the Hough functions depend is 

q = 2Sl*/a km , 

which justifies writing H\ m and e\i m rather than H\ m and e\ lm . 



(A12) 



(A13) 



We again represent the linear response of a star, as in Appendix Al by a vector y(r, 0, <f>, i) whose components are the 
various oscillation variables, and which can be expressed in the inertial frame as: 



V = 



M 
Mi 

Ah 
Mi 



exp(— ikQ olb t) ^2 

k I 

exp(— ikQ olb t) ^ 



" W lm e^xL(e) £ H k Xm (n)y k xlm {r) 

TTl A 



(A14) 



The expansion of H Xm back into associated Legendre functions in the second line of eq. (A14l is useful since disk integrals 
are convenient to perform over spherical harmonics (§ 3.31. 

Following | Unno et al.| ( |1989[ ), we choose the components of y as 



2/1 



6- 



2/2 



Sp 
pgr' 



As , AL 

— , and ye = — , 



(A15) 



where we have omitted the variables corresponding to the perturbed gravitational potential, yz and j/4. Equation (A14l 
together with determination of the radial displacement £ r /r = y 1 and the Lagrangian flux perturbation AF/F = ye — 2yi at 
the photosphere then enables use of the formalism from § |3.3| to compute the flux perturbation as seen by an observer. 

Next, we present the differential equations which determine a particular component y X i m { r ) °f t ne full response in 
radiative zones. These equations are nearly identical to those in the appendix of Pfahl et al. ( 2008 1, but with 1(1 + 1) replaced 
by A and with certain terms set to zero as per the traditional approximation. In practice these terms can be left in, since they 
are nearly zero for situations where the traditional approximation is valid; this is then a smooth way of transitioning among 
different regimes. Omitting (Xlmk) indices and denoting U = U\ lm and w = akm, the equations are 



dyi 
dlnr 

dy2 
dlnr 

dys 
dlnr 



= J/i 



= 2/i 



2/i 



3 +2/2 



Hu 



gr 



uj 2 -N' 

g/r 

Vad ( r? 



As 



gr 



uj 2 r^ 



2/5 Ps + T^2 U 



+ 2/2 1 - T) 



N 2 



ysPs 



g/r 



g/r 

+ 4(V- V ad ) + c 2 



1 dU 
g dr 
r 



+ 2/2 



H„ 



(V a 



+ ys ir V(4 : - k s 

ti. 



r „ r 
?/6— V + 



dye 
dlnr 



2/2 



(¥) 



+ 2/5 [iui 



4nr 



''P£pT\ 



2/67 ' 



dU/dr 



+ (V ad - V)- 



A 



u 



u. 



(A16) 
(A17) 

(A18) 

(A19) 



where r) = 47rr 3 p/M r , 7 = 47rr 3 pe/L r , c 2 = (r/ff p )V(K ad - 4V a 
scale height, and e is the specific energy generation rate. 



V a d(rfln V a d/rflnr + r/H p ), H p = pg/p is the pressure 
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We need four boundary conditions for our four variables. Our first three are 

= £r(0) evanescence in convective core, (A20) 

= As(0) adiabaticity/evanescence in core, (A21) 

= ~E7^T ~ 4 ~PFT^T blackbody at the stellar surface, (A22) 

F(R) T(R) 

where AT/T can be cast in terms of yi, 3/2, and j/5 using standard thermodynamic derivative identities. 

A final surface boundary condition that allows for traveling and/or standing waves can be generated by imposing energetic 
constraints at the surface. This is detailed in |Unno et al.| ( fl989[ l pp. 163 - 167 for adiabatic oscillations. To generalize the 



boundary condition to include nonadiabaticity, rotation, and inhomogeneous tidal forcing, we write equations (A16l - (A19l 

as 

^ r =My + b , (A23) 

where M and b are treated as constant near the stellar photosphere. The constant solution is y — —M~ 1 b; defining 
z = y — y , the homogeneous solutions for z can be computed by diagonalizing M. In the evanescent case, we eliminate the 
solution for y with outwardly increasing energy density. Alternatively, in the traveling wave case, we eliminate the inward- 
propagating wave. The final boundary condition is then implemented by setting the amplitude of the eliminated homogeneous 
solution to zero, and solving for a relationship between the original fluid variables implied by this statement. 



APPENDIX B: ANALYTIC MODEL OF ELLIPSOIDAL VARIATION 

As discussed in § |6.1[ our simplified model of ellipsoidal variation reproduces the much more sophisticated simulation code 



employed by Welsh et al. (20111 to model KOI-54; here we discuss the details of our analytic methods, which can easily be 



applied to model other systems. 



Bl Irradiation 

The following is our simple analytical model of the insolation component of the KOI-54's ellipsoidal variation. We focus our 
analysis on the primary, since extending our results to the secondary is trivial. Our main assumption is that all radiation from 
the secondary incident upon the primary is immediately reprocessed at the primary's photosphere and emitted isotropically 
(i.e., absorption, thermalization, and reemission). This assumption is well justified for KOI-54, since its two component stars 
are of very similar spectral type. The method below might need to be modified if the components of a binary system had 
significantly different spectral types, because then some of the incident radiation might instead be scattered. 
The incident flux on the primary, using the conventions and definitions introduced in § |3.1| is 



where Z is the ramp function, defined by 



F2 - = 4^^-^ (B1) 



*(»)= < ° X< " ■ (B2) 
x x > 



We can expand Z[f ■ D) in spherical harmonics as 

Z(r-D) = J2 Zl ™ Y >™( 6 >^ e ~ lmm > ( B3 ) 



where Zi m can be evaluated to 



/ 211 + 1 (l-m)!\ 1/2 /cos(m7r/2)\ f 1 , r- , 



(B4) 



with cos(m7r/2)/(l — m 2 ) — > 7r/2 for m — ±1. 

Next, taking the reemission as isotropic, the reemitted intensity will be 



/emit = — • (B5) 



7T 



Using this together with our expansion of Z(f-D) as well as results from § 3.3 we can evaluate the observed flux perturbation 



5 2 \ v 1 



T=^)[T--imt)T, E hZ lm Y lm {e o ,<t>o)e~ im ^\ (B6) 



1=0 m=-l 

„2 



where Ji = Li/ins is the unperturbed observed flux, s is the distance to the observer, the bandpass correction coefficient 
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/3(T) is defined in eq. ( |19[ ), the disk-integral factor bi is denned in eq. (231, several values of bi using Eddington limb darkening 
are given in Table |3j and other variables are defined in § |3.1| Since bi declines rapidly with increasing I, it is acceptable to 
include only the first few terms of the sum in eq. (B6l. We have neglected limb darkening, so it is formally necessary to use 
a flat limb darkening law in calculating bi (h(fj,) = 2; § 3.3 1. However, we found this to be a very modest effect. 

The binary separation D and the true anomaly / can be obtained as functions of time in various ways, e.g., by expanding 
with the Hansen coefficients employed earlier (eq. 11 or A6l, or by using 

D= a ^- e2 \ (B7) 



1 + e cos / 



together with numerical inversion of 



fiorbi = 2 arctan 



(l-e)tan(//2) 



! sin/ 



1 + e cos / 



-n < f < tv. 



(B8) 



+ 7T. 



The observed flux perturbation from the secondary is obtained from eq. ( B6 1 by switching If) 2 and sending (f> - 

Using the fact that bo = 1 and Zoo = \A"/2, it can be readily verified that the total reflected power, i.e. s 2 Ji times 
( B6 1 integrated over all observer angles (0 O , <j> ), is equal to Lz^R^/AnD 2 ). This is just the secondary's luminosity times the 
fraction of the secondary's full solid angle occupied by the primary, which is the total amount of the secondary's radiation 
incident on the primary. 



B2 Equilibrium tide 

We invoke the Cowling approximation (well satisfied for surface values of perturbation variables), and use the analytic 
equilibrium tide solution, where the radial displacement at star l's surface becomes 

U(Ri,t) 



(B9) 



19891 



and U is the tidal potential. Using the expansion in eq. (A2|, £ r ,im(t)/Ri from eq. (201 becomes (Goldreich fc Nicholson) 



jr,l m {t) = M2 / Rl 

III Mi \D(t) 



Wi m e" 



(BIO) 



We invoke von Zeipel's theorem ( von Zeipel 1924 Pfahl et al. 2008 1 to determine the corresponding surface emitted flux 
perturbation: 



AF lm {t) 



-(1 + 2) 



We can then explicitly evaluate the observed flux variation using the formalism from § |3.3| 
Ji Mi 



' )J U2 zZ(ms) E {[2-l3(Ti)(l + 2))b l -c l }w lm Y lm (e o ,<p a ) 



1=2 v v ' ' m=-l 



-imf(t) 



(Bll) 



(B12) 



where the bandpass correction coefficient j3(T) is defined in eq. (19 1, Wi m is defined in eq. (A3 1, the disk-integral factors bi 
and c; are defined in equations (231 and (24 1, several values of bi and c\ using Eddington limb darkening are given in Table [3j 
and other variables are defined in § |3.1| Due to the strong dependence on I, it is typically acceptable to include only the first 
term of the sum in eq. (B6l. Computation of the binary separation D(t) and true anomaly f(t) is discussed in Appendix [Bl] 
The observed flux perturbation from the secondary is obtained from eq. (B12| by switching 1 o 2 and sending (j> — > <j) + n. 

We note that although the analytic equilibrium tide solution for the radial displacement £ r is a good approximation at 
the stellar surface regardless of stellar parameters, the presence of a significant surface convection zone in a solar-type star 
proscribes the use of eq. (Bll I; Pfahl et al. ( 2008 1 gives the appropriate replacement in their eq. (37). Moreover, we note that 
eq. (JbTTJ) may also be invalid for slowly rotating stars in eccentric orbits; see § |6.2| 



APPENDIX C: TIDAL ORBITAL EVOLUTION 

CI Eigenmode expansion of tidal torque and energy deposition rate 

Assuming alignment of rotational and orbital angular momenta, the tidal torque r produced by star 2 on star 1 must have 



only a z component, where z points along the orbital angular momentum. We can evaluate it as follows (Kumar & Quataert 



19981. First, 



IX r ' d ^) iV "L {i " ) 

4>- [-(po + 8p)VU]r sin9 dV. 



dF 

dV 



dV 



(CI) 
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The term involving the background density po vanishes; expanding both the tidal potential U and the Eulerian density 
perturbation 5p in spherical harmonics with expansion coefficients Uim(r,t) and 8pi m (r,t), we have 



oo 1 - Bl 

T { t ) = *E E 711 / bPl™{r 
1 = 2 m = -l 



t)Ui m (r, t)r 2 dr. 



(C2) 



Further invoking the expansions from equations (|5| and (A4|, as well as the definitions in § 3.2 we arrive with 



i (£iQn;W; m ) 2 ~fc ~fe' A i(fc'-fc)n orb t 



m - -" XLXf m A wimfc e^ -^w (C3) 

Lastly, averaging over a complete orbital period and rearranging the sums, we derive our final expression for the secular tidal 
torque: 



(r} = 8 



( GMt 
V Ri 



M2 
Mi 



E 



21+2 1 



(C4) 



The torque depends on the rotation rate fi* only through the Doppler-shifted frequency akm = fc^orb — mfi* , since we have 
neglected rotational modification of the eigenmodes (§ 3.4 1. Fig. [4] shows plots of this torque evaluated for KOF54. 

Note that a particular term of this sum is positive if and only if matm ~ m(fcf2 or b — mfl.) > 0, which reduces to 
(k/rn)tt OI b > O*. This is known as being prograde, since it is equivalent to the condition that a mode's angular structure, in 
the corotating frame, rotate in the same sense as the stellar spin; conversely, retrograde waves with (fc/m)fi or b < n» cause 
negative torques. 

Using similar techniques to those given above, an equivalent expansion of the secular tidal energy deposition rate into 
the star (including mechanical rotational energy) can be derived: 

21+2 1 



(E) = 8fi rb 



( GMl 



Ma 
Mi 



E»i.E«t.WE(g)( 

rr> — — / Ir — 11 n x ' N 



1 = 2 v ' m=—l k=0 

The only difference between equations |C4l and (|C5h is switching m o fcfiorb- 



2 _ a 2 \ 2 _i_ A-,2 2 
ril km' 1 ' rtl km 



(C5) 



C2 Nonresonant pseudosynchronization 

A pseudosynchronous frequency Sl ps is defined as a rotation rate that produces no average tidal torque on the star throughout 
a sufficiently long time interval, which here we take to be a complete orbital period (§[H|. I.e., 

(T>(Op.)=0. (C6) 



Here we will show that our expansion from CI reproduces the value of fi ps derived in Hut ( 1981 1, which we denote ttpl, in the 
equilibrium tide limit. We will in particular show that Hut's result is independent of assumptions about eigenmode damping 
rates. 

Proceeding, we take the nonresonant (equilibrium tide) limit of eq. (C4|. This is obtained by retaining only the first term 
in the Taylor series expansion in crfc m /w„; of the last factor in parentheses from eq. (|C4[), and yields 



(Tnr) = 8 



(GMf \ (M 2 \ 2 ^(R^ 2l+2 



V Ri 



Mi 



E 



E 



Qnl \ ( IrU 

Em) U 2 , 



(C7) 



note that sums over k and m become decoupled from the sum over n. Setting (r nr )(f2ps) = and retaining only I — 2, we 
have 



0= J2 A 2 fc 2 (e) 2 (fcfi rb-2^). 

k = — oo 

We need two identities to evaluate this further. First, starting with the definition of the Hansen coefficients, 



'a\t +1 -imf 



- E * 



k — ifcr2 or ^,t 
lm ^ j 



(C8) 



(C9) 



we can differentiate with respect to t, then multiply by the complex conjugate of ( C9 1 and average over a complete period to 
derive 



00 



2 m 
~ 2^ 



1 + e cos / 
1-e 2 



2i+2 



df. 



(CIO) 



Specializing to I = 2, 



E *(* 



rk 
L 2lTl 



5e 6 + 90e 4 + 120e 2 + 16 



16(1 -e 2 ) 6 



(Cll) 
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The second identity needed, 



3e 4 + 24e 2 + 8 

g ( 1 _ e 2)3/2 ' 



(C12) 



can be derived similarly. 



Substituting equations (Clll and ( C12 1 into |C8|, we have that 



nr _ „ 1 + (15/2)e 2 + (45/8)e 4 + (5/16)e 
ps ~ ° rb ' [1 + 3e 2 + (3/8)e 4 ] (1 - e 2 ) 3 / 2 



(C13) 



this is precisely eq. (42) from Hut ( 1981 1 
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